Run Summary

import pandas as pd
print(df.head())
print(df.columns)
print(df.shape)
print(df.describe())

print(df['star1_x'].mean(), df['star1_y'].mean(), df['star2_x'].mean(), df['star2_y'].mean())

print(df[['star2_x','star2_y','star2_z']].head())

print(df.iloc[0])

import numpy as np
# Compute separation between stars for all time
r_sep = np.sqrt((df['star1_x'] - df['star2_x'])**2 + (df['star1_y'] - df['star2_y'])**2 + (df['star1_z'] - df['star2_z'])**2)
print(r_sep.mean(), r_sep.std(), r_sep.std()/r_sep.mean())

# Let's compute separation and peak times
import numpy as np
sep = np.sqrt((df['star1_x'] - df['star2_x'])**2 + (df['star1_y'] - df['star2_y'])**2 + (df['star1_z'] - df['star2_z'])**2).values
times = df['time'].values
# find local maxima: indices where sep[i]>sep[i-1] and sep[i]>sep[i+1]
peaks = np.where((sep[1:-1] > sep[:-2]) & (sep[1:-1] > sep[2:]))[0] + 1
# get first few peaks
peak_times = times[peaks]
peak_seps = sep[peaks]
# print first few
print(len(peaks))
print('First few peak times & separations:')
for i in range(5): print(f"peak {i}: time = {peak_times[i]}, sep = {peak_seps[i]}")
# compute time differences between successive peaks
dts = peak_times[1:] - peak_times[:-1]
# print some dt
print('First few dt: ', dts[:5])
print('Mean dt:', np.mean(dts), 'Median dt:', np.median(dts), 'Std dt:', np.std(dts))
# use median dt as period
dt_median = np.median(dts)
print('Period ~', dt_median, 's')

print(sep.min(), sep.max())

print([v for v in globals() if not v.startswith('_')])

import numpy as np
# Given values
g = 6.67430e-11
r_min = sep.min()
r_max = sep.max()
a = (r_min + r_max) / 2
dt_median = dt_median # from previous code
# Compute total mass using Kepler's third law
total_mass = 4 * np.pi**2 * a**3 / (g * dt_median**2)
# Compute total energy E = - G M1 M2 / (2a), with M1=M2=total_mass/2
# So E = -G(total_mass/2 * total_mass/2)/(2a) = -G total_mass^2/(8a)
E_total = - g * total_mass**2 / (8 * a)
print("r_min =", r_min)
print("r_max =", r_max)
print("a =", a)
print("Period P =", dt_median)
print("Total mass (kg) =", total_mass)
print("Total mass (sol mass) =", total_mass / 1.98847e30)
print("Total energy (J) =", E_total)

# Compute mass ratio M1/M2 from ratio of velocity magnitudes v2/v1
import numpy as np
# positions
r1 = df[['star1_x','star1_y','star1_z']].values
r2 = df[['star2_x','star2_y','star2_z']].values
# time
t = df['time'].values
# compute displacement vectors for central difference
dr1 = r1[2:] - r1[:-2]
dr2 = r2[2:] - r2[:-2]
# times dt
dt = t[2:] - t[:-2]
# velocity magnitude approximations
dv1 = dr1/ dt[:,None]
dv2 = dr2/ dt[:,None]
speed1 = np.linalg.norm(dv1, axis=1)
speed2 = np.linalg.norm(dv2, axis=1)
# ratio of speeds
ratio = speed2 / speed1
# filter out unrealistic values, e.g. speed1 or speed2 extremely small or ratio inf/nan
valid = np.isfinite(ratio) & (speed1>0) & (speed2>0)
ratio_valid = ratio[valid]
# take median ratio and mean ratio
median_ratio = np.median(ratio_valid)
mean_ratio = np.mean(ratio_valid)
print('median M1/M2:', median_ratio)
print('mean M1/M2:', mean_ratio)
print('some stats:', np.percentile(ratio_valid, [1,5,25,50,75,95,99]))
print('number of valid points:', len(ratio_valid))

# Compute energies using instantaneous velocities and computed masses
import numpy as np
# Given gravitational constant
G = 6.67430e-11
# Period P and semi-major axis a (recomputed here for clarity)
...
# But we'll define P and a from previous code:
P = dt_median
a = (r_min + r_max) / 2
# Total mass from third law
M_total = 4 * np.pi**2 * a**3 / (G * P**2)
# Mass ratio q = median of speed2/speed1
q = median_ratio
# Compute individual masses
M1 = q / (1 + q) * M_total
M2 = 1 / (1 + q) * M_total
print("M_total, M1, M2 (kg): ", M_total, M1, M2)
# Compute velocities via central difference dr/dt from r arrays
dt = df['time'].values
r1 = df[['star1_x','star1_y','star1_z']].values
r2 = df[['star2_x','star2_y','star2_z']].values
# central diff
dr1 = r1[2:] - r1[:-2]
dr2 = r2[2:] - r2[:-2]
dt_c = dt[2:] - dt[:-2]
vel1 = dr1 / dt_c[:, np.newaxis]
vel2 = dr2 / dt_c[:, np.newaxis]
# speeds
speed1 = np.linalg.norm(vel1, axis=1)
speed2 = np.linalg.norm(vel2, axis=1)
# separation for central points r_sep_c
dr_rel = r1[1:-1] - r2[1:-1]
r_sep_c = np.linalg.norm(dr_rel, axis=1)
# energies
K1 = 0.5 * M1 * speed1**2
K2 = 0.5 * M2 * speed2**2
U = - G * M1 * M2 / r_sep_c
E = K1 + K2 + U
print('mean E:', np.mean(E), 'std E:', np.std(E))
print('E based on analytic formula:', -G*M1*M2/(2*a))

{'answer': -1.8785837068999796e+37}

import pandas as pd
print(df.columns)
print(df.head())

import numpy as np
# compute r1, r2 and r_rel magnitudes
df['r1'] = np.sqrt(df['star1_x']**2 + df['star1_y']**2)
df['r2'] = np.sqrt(df['star2_x']**2 + df['star2_y']**2)
df['r_rel'] = np.sqrt((df['star1_x']-df['star2_x'])**2 + (df['star1_y']-df['star2_y'])**2)
print('r1_stats:', df['r1'].min(), df['r1'].max())
print('r2_stats:', df['r2'].min(), df['r2'].max())
print('r_rel_stats:', df['r_rel'].min(), df['r_rel'].max())

print(df['time'].min(), df['time'].max())

print(df['r1'].mean(), df['r2'].mean())

# Check if r1 + r2 equals r_rel
import numpy as np
diff = df['r1'] + df['r2'] - df['r_rel']
print(diff.min(), diff.max(), np.mean(diff))

import numpy as np
import pandas as pd
# calculate r_rel
df['x_rel'] = df['star1_x'] - df['star2_x']
df['y_rel'] = df['star1_y'] - df['star2_y']
df['r_rel'] = np.sqrt(df['x_rel']**2 + df['y_rel']**2)
# find local minima and maxima in r_rel
good = np.arange(1, len(df)-1)
mins = []
maxs = []
for i in good:
    if df['r_rel'].iloc[i-1] > df['r_rel'].iloc[i] < df['r_rel'].iloc[i+1]:
        mins.append((df['time'].iloc[i], df['r_rel'].iloc[i], i))
    if df['r_rel'].iloc[i-1] < df['r_rel'].iloc[i] > df['r_rel'].iloc[i+1]:
        maxs.append((df['time'].iloc[i], df['r_rel'].iloc[i], i))
# convert to arrays
mins = np.array(mins, dtype=object)
maxs = np.array(maxs, dtype=object)
# print some
print('Found minima:', len(mins), ' maxima:', len(maxs))
# get few
print('First 5 mins:', mins[:5])
print('First 5 maxs:', maxs[:5])
# compute period differences from minima
if len(mins)>1:
    diff_times = np.diff(mins[:,0].astype(float))
    P = np.median(diff_times)
    print('Estimated period P (s):', P, ' = days', P/86400)
# compute Rmin and Rmax median
if len(mins)>0 and len(maxs)>0:
    Rmin = np.median(mins[:,1].astype(float))
    Rmax = np.median(maxs[:,1].astype(float))
    print('Estimated r_min (m):', Rmin)
    print('Estimated r_max (m):', Rmax)
    a = (Rmin + Rmax) / 2
    print('Semi-major a (m):', a)
    e = (Rmax - Rmin) / (Rmax + Rmin)
    print('Eccentricity e:', e)

# compute ratio r1/r2 and check variation
ratio = df['r1']/df['r2']
print(ratio.mean(), ratio.std(), ratio.min(), ratio.max())

# constants
i = df
import numpy as np
G = 6.67430e-11
P = 430968429.988678
a = 796777638309.3845
total_mass = 4*np.pi**2 * a**3 / (G * P**2)
print('M_total:', total_mass, 'kg')

print(df.shape)

import numpy as np
# Downsample data to speed up SSE scans
downsample_rate = 100
i_idxs = np.arange(0, len(df), downsample_rate)
t = df['time'].values[i_idxs]
x1 = df['star1_x'].values[i_idxs]
y1 = df['star1_y'].values[i_idxs]
x2 = df['star2_x'].values[i_idxs]
y2 = df['star2_y'].values[i_idxs]
# precompute variables (for performance), but we will use polyfit in loop
def sse_for_k(k_val):
    # compute COM coords
    COM_x = (k_val * x1 + x2) / (k_val + 1)
    COM_y = (k_val * y1 + y2) / (k_val + 1)
    # fit line COM_x = ax*t + bx
    a_x, b_x = np.polyfit(t, COM_x, 1)
    # fit line COM_y = ay*t + by
    a_y, b_y = np.polyfit(t, COM_y, 1)
    # compute SSE
    sse_x = np.sum((COM_x - (a_x * t + b_x))**2)
    sse_y = np.sum((COM_y - (a_y * t + b_y))**2)
    return sse_x + sse_y
# coarse grid scan
ks = np.logspace(-2, 2, 200) # from 0.01 to 100
sses = np.array([sse_for_k(k) for k in ks])
# find best k
i_best = np.argmin(sses)
k_best_coarse = ks[i_best]
print('Coarse best k:', k_best_coarse, 'SSE:', sses[i_best])
# refine around k_best_coarse using local optimization
from scipy.optimize import minimize_scalar
# define objective for minimize_scalar (only needs positive k)
res = minimize_scalar(lambda k: sse_for_k(k), bracket=[ks[max(i_best-5,0)], ks[min(i_best+5,len(ks)-1)]], bounds=(1e-3,1e3), method='bounded', options={'xatol':1e-6})
k_best = res.x
print('Refined best k:', k_best, 'SSE:', res.fun)

# calculate total energy
g = 6.67430e-11
M_total = 1.6109219860276977e30
k = 0.28571427624283185
a = 796777638309.3845
m1 = k/(1+k) * M_total
m2 = 1/(1+k) * M_total
E_total = -g*m1*m2/(2*a)
print("m1 =", m1, "kg")
print("m2 =", m2, "kg")
print("Total energy E =", E_total, "J")

# Confirm k_best with downsample_rate=50
import numpy as np, pandas as pd
df2 = df.copy()
downsample_rate = 50
i_idxs = np.arange(0, len(df2), downsample_rate)
t = df2['time'].values[i_idxs]
x1 = df2['star1_x'].values[i_idxs]
y1 = df2['star1_y'].values[i_idxs]
x2 = df2['star2_x'].values[i_idxs]
y2 = df2['star2_y'].values[i_idxs]
def sse_for_k(k_val):
    COM_x = (k_val * x1 + x2) / (k_val + 1)
    COM_y = (k_val * y1 + y2) / (k_val + 1)
    # linear fit to COM coords
a_x, b_x = np.polyfit(t, COM_x, 1)
a_y, b_y = np.polyfit(t, COM_y, 1)
def sse_fun(k):
    COM_x = (k * x1 + x2) / (k + 1)
    COM_y = (k * y1 + y2) / (k + 1)
    sse_x = np.sum((COM_x - (a_x * t + b_x))**2)
    sse_y = np.sum((COM_y - (a_y * t + b_y))**2)
    return sse_x + sse_y
from scipy.optimize import minimize_scalar
res = minimize_scalar(sse_fun, bounds=(1e-3, 1e3), method='bounded', options = {'xatol':1e-6})
print('k_best50:', res.x)

{'answer': -1.8785839928397858e+37}

import numpy as np
# let's inspect the DataFrame structure
df.head()

print(df.head())

print(df.columns)

print({col: df.iloc[0][col] for col in df.columns})

print(df.time.min(), df.time.max())

print(df[['time','r_rel']].tail())
print(df.r_rel.min(), df.r_rel.max())

print("star1_x: min, max", df.star1_x.min(), df.star1_x.max())
print("star1_y: min, max", df.star1_y.min(), df.star1_y.max())
print("star2_x: min, max", df.star2_x.min(), df.star2_x.max())
print("star2_y: min, max", df.star2_y.min(), df.star2_y.max())

import numpy as np
r_rel = df['r_rel'].values
times = df['time'].values
# find minima
good = (r_rel[1:-1] < r_rel[:-2]) & (r_rel[1:-1] < r_rel[2:])
minima = np.where(good)[0] + 1
times_min = times[minima]
diffs = np.diff(times_min)
print("Number of minima: ", len(times_min))
print("Mean period estimate (s):", np.mean(diffs), "std:", np.std(diffs))
print("First few minima times:", times_min[:10])
print("First few period diffs:", diffs[:10])

import numpy as np
G = 6.67430e-11
# DataFrame df is loaded
r_rel = df['r_rel'].values
times = df['time'].values
# find local minima indices for r_rel
good = (r_rel[1:-1] < r_rel[:-2]) & (r_rel[1:-1] < r_rel[2:])
minima_idx = np.where(good)[0] + 1
# compute period from differences of minima times
P = np.mean(np.diff(times[minima_idx]))
# compute semimajor axis a
r_min = np.min(r_rel)
r_max = np.max(r_rel)
a = (r_min + r_max) / 2
# total mass from Kepler's third law
total_mass = 4 * np.pi**2 * a**3 / (G * P**2)
# mass ratio from barycenter distances: m1/m2 = r2/r1 constant	ratio_m1_to_m2 = np.mean(df['r2'] / df['r1'])
# solve for individual masses
# m1/m2 = ratio, so m1 = ratio * m2, so total_mass = ratio*m2 + m2 = m2*(1+ratio)
# => m2 = total_mass/(1+ratio), m1 = total_mass - m2
m2 = total_mass / (1 + ratio_m1_to_m2)
m1 = total_mass - m2
# total energy E = -G m1 m2 / (2a)
E_total = -G * m1 * m2 / (2 * a)
print("Computed values:")
print("P (s):", P)
print("r_min (m):", r_min)
print("r_max (m):", r_max)
print("a (m):", a)
print("total_mass (kg):", total_mass)
print("m1/m2 ratio:", ratio_m1_to_m2)
print("m1 (kg):", m1)
print("m2 (kg):", m2)
print("E_total (J):", E_total)

import numpy as np;
G=6.67430e-11
r_rel = df['r_rel'].values
times = df['time'].values
good = (r_rel[1:-1] < r_rel[:-2]) & (r_rel[1:-1] < r_rel[2:])
minima_idx = np.where(good)[0] + 1
P_vals = np.diff(times[minima_idx])
P = np.mean(P_vals)
P_std = np.std(P_vals)
r_min = np.min(r_rel)
r_max = np.max(r_rel)
a = (r_min + r_max) / 2
total_mass = 4 * np.pi**2 * a**3 / (G * P**2)
ratio = np.mean(df['r2'] / df['r1'])
m2 = total_mass / (1 + ratio)
m1 = total_mass - m2
E_total = -G * m1 * m2 / (2 * a)
print("P_vals first few:", P_vals[:5])
print("P mean, P std:", P, P_std)
print("r_min, r_max, a:", r_min, r_max, a)
print("total_mass:", total_mass)
print("ratio:", ratio)
print("m1, m2:", m1, m2)
print("E_total:", E_total)

# Let's compute instantaneous energies and verify
import numpy as np
G = 6.67430e-11
# Use previously computed m1, m2, a
a = 7.967776111190277e11
m1 = 7.268624999400248e29
m2 = 8.841692835765647e29
# positions
t = df['time'].values
x1 = df['star1_x'].values
y1 = df['star1_y'].values
x2 = df['star2_x'].values
y2 = df['star2_y'].values
r_rel = df['r_rel'].values
# velocities via central difference
tol = 10  # sample size to check some points
ii = np.linspace(10000, len(t)-10000, num=10, dtype=int)
energies = []
for i in ii:
    # velocity star1
dt1 = t[i+1] - t[i-1]
v1x = (x1[i+1] - x1[i-1]) / dt1
v1y = (y1[i+1] - y1[i-1]) / dt1
v1_sq = v1x**2 + v1y**2
# velocity star2
dt2 = t[i+1] - t[i-1]
v2x = (x2[i+1] - x2[i-1]) / dt2
v2y = (y2[i+1] - y2[i-1]) / dt2
v2_sq = v2x**2 + v2y**2
# KE
e1 = 0.5 * m1 * v1_sq
n2 = 0.5 * m2 * v2_sq
# potential energy
U = -G * m1 * m2 / r_rel[i]
energies.append(e1 + n2 + U)
energies = np.array(energies)
print("Energies sample:", energies)
print("Mean energy:", energies.mean(), "std:", energies.std())

print("r1 min,max", df['r1'].min(), df['r1'].max())
print("r2 min,max", df['r2'].min(), df['r2'].max())

import numpy as np
# Constants
g = 6.67430e-11
# Masses and semi-major axis from earlier
m1 = 7.268624999400248e29
m2 = 8.841692835765647e29
# Data
t = df['time'].values
x1 = df['star1_x'].values
y1 = df['star1_y'].values
x2 = df['star2_x'].values
y2 = df['star2_y'].values
r_rel = df['r_rel'].values
# choose sample indices avoiding edges
i_samples = np.linspace(10000, len(t)-10000, num=20, dtype=int)
energies = []
for i in i_samples:
    # central difference for velocities
    dt = t[i+1] - t[i-1]
    v1x = (x1[i+1] - x1[i-1]) / dt
    v1y = (y1[i+1] - y1[i-1]) / dt
    v2x = (x2[i+1] - x2[i-1]) / dt
    v2y = (y2[i+1] - y2[i-1]) / dt
    v1_sq = v1x**2 + v1y**2
    v2_sq = v2x**2 + v2y**2
    K = 0.5 * m1 * v1_sq + 0.5 * m2 * v2_sq
    U = -g * m1 * m2 / r_rel[i]
    energies.append(K + U)
energies = np.array(energies)
print("Sample energies:", energies)
print("Mean energy:", energies.mean())
print("Std energy:", energies.std())
print("Expected formula energy:", -g * m1 * m2 / (2 * ((r_rel.min() + r_rel.max())/2)))

import numpy as np;
G=6.67430e-11
r_rel = df['r_rel'].values
t = df['time'].values
# find minima times for r_rel
good = (r_rel[1:-1] < r_rel[:-2]) & (r_rel[1:-1] < r_rel[2:])
minima_idx = np.where(good)[0] + 1
times_min = t[minima_idx]
# refine period
dt_period = (times_min[-1] - times_min[0]) / (len(times_min) - 1)
# a from separation extremes
r_min = np.min(r_rel)
r_max = np.max(r_rel)
a = (r_min + r_max) / 2
# total mass
total_mass = 4 * np.pi**2 * a**3 / (G * dt_period**2)
ratio = np.mean(df['r2'] / df['r1'])
m2 = total_mass / (1 + ratio)
m1 = total_mass - m2
E = -G * m1 * m2 / (2 * a)
print("P0 period refined:", dt_period)
print("a:", a)
print("total_mass:", total_mass)
print("ratio m1/m2:", ratio)
print("m1, m2:", m1, m2)
print("E_total refined:", E)

# Step 1–3: find period P, semimajor axis a, total mass M
def compute_orbital_parameters(df):
    import numpy as np
    from scipy.signal import argrelextrema
    G = 6.67430e-11
    t = df['time'].values
    r_rel = df['r_rel'].values
    # find all local minima in r_rel
    # we can use a simple comparison
    minima = (r_rel[1:-1] < r_rel[:-2]) & (r_rel[1:-1] < r_rel[2:])
    idx_min = np.where(minima)[0] + 1
    # require at least 2 minima
    if len(idx_min) < 2:
        raise ValueError("Not enough periapses found")
    t_min = t[idx_min]
    # orbital period estimate
    P_vals = np.diff(t_min)
    P = np.mean(P_vals)
    # a from rmin, rmax
    rmin = np.min(r_rel[idx_min])
    # for rmax use local maxima too
    maxima = (r_rel[1:-1] > r_rel[:-2]) & (r_rel[1:-1] > r_rel[2:])
    idx_max = np.where(maxima)[0] + 1
    rmax = np.max(r_rel[idx_max])
    a = 0.5*(rmin + rmax)
    # total mass from Kepler
    M = 4*np.pi**2 * a**3 / (G * P**2)
    return P, a, M, idx_min, idx_max

P, a, M, idx_min, idx_max = compute_orbital_parameters(df)
print(f"P = {P:.3e} s, a = {a:.3e} m, M = {M:.3e} kg")

import numpy as np
G = 6.67430e-11
df2 = df.copy()
# compute local minima and maxima
r_rel = df2['r_rel'].values
t = df2['time'].values
# local minima indices
minima = (r_rel[1:-1] < r_rel[:-2]) & (r_rel[1:-1] < r_rel[2:])
idx_min = np.where(minima)[0] + 1
# local maxima indices
maxima = (r_rel[1:-1] > r_rel[:-2]) & (r_rel[1:-1] > r_rel[2:])
idx_max = np.where(maxima)[0] + 1
r_min_vals = r_rel[idx_min]
r_max_vals = r_rel[idx_max]
t_min_vals = t[idx_min]
# refine period P
P_vals = np.diff(t_min_vals)
P = np.mean(P_vals)
P_std = np.std(P_vals)
# a
a_local = 0.5 * (np.mean(r_min_vals) + np.mean(r_max_vals))
# total mass
M_local = 4 * np.pi**2 * a_local**3 / (G * P**2)
# ratio m1/m2
ratio = np.mean(df2['r2'].values / df2['r1'].values)
m2 = M_local / (1 + ratio)
m1 = M_local - m2
E = -G * m1 * m2 / (2 * a_local)
print("Local minima (r) average:", np.mean(r_min_vals), "std:", np.std(r_min_vals))
print("Local maxima (r) average:", np.mean(r_max_vals), "std:", np.std(r_max_vals))
print("P: mean", P, "std", P_std)
print("a_local:", a_local)
print("M_local:", M_local)
print("ratio m1/m2:", ratio)
print("m1 (kg):", m1)
print("m2 (kg):", m2)
print("E_total (J):", E)

{'answer': -2.691698195548385e+37}

import pandas as pd
# show columns
df.columns
df.head()

print(df.columns)
print(df.head())

print(df.iloc[0])

print(df.shape)

print(globals().keys())

import numpy as np
# positions
df_np = df
r1 = df_np[['star1_x','star1_y','star1_z']].values
r2 = df_np[['star2_x','star2_y','star2_z']].values
# times
t = df_np['time'].values
# relative vector
r_vec = r2 - r1
r_mag = np.linalg.norm(r_vec, axis=1)
# angles and period
phi = np.arctan2(r_vec[:,1], r_vec[:,0])
phi_un = np.unwrap(phi)
delta_phi = phi_un[-1] - phi_un[0]
cycles = delta_phi/(2*np.pi)
T = (t[-1] - t[0]) / cycles
# separation extremes
r_max = r_mag.max()
r_min = r_mag.min()
a = (r_max + r_min)/2
# total mass
G = 6.67430e-11
M_total = 4*np.pi**2 * a**3 / (G * T**2)
# velocities using central difference
# dt assume constant;
but use vector dt_i average
# Actually t spacing constant
# Use dt = t[i+1] - t[i]
dt = t[1] - t[0]
v1 = (r1[2:] - r1[:-2])/(2*dt)
v2 = (r2[2:] - r2[:-2])/(2*dt)
v1_mag = np.linalg.norm(v1, axis=1)
v2_mag = np.linalg.norm(v2, axis=1)
ratio = v2_mag / v1_mag
mass_ratio = np.median(ratio)
m1 = M_total * mass_ratio/(1 + mass_ratio)
m2 = M_total - m1
E_total = -G * m1 * m2 / (2 * a)
# print results
print("Total time span (s):", t[-1] - t[0])
print("Number of cycles:", cycles)
print("Orbital period (s):", T)
print("Min separation (m):", r_min)
print("Max separation (m):", r_max)
print("Semi-major axis (m):", a)
print("Total mass (kg):", M_total)
print("Mass ratio m1/m2:", mass_ratio)
print("Mass1 (kg):", m1)
print("Mass2 (kg):", m2)
print("Total energy (J):", E_total)

# Let's refine code: compute everything carefully.
import numpy as np
# get positions and times
t = df['time'].values
r1 = df[['star1_x','star1_y','star1_z']].values
r2 = df[['star2_x','star2_y','star2_z']].values

# relative vector r12 = r2 - r1
def compute_period_and_a(t, r1, r2):
    r_vec = r2 - r1
    # magnitude
    r_mag = np.linalg.norm(r_vec, axis=1)
    # extremums: r_min, r_max
    r_min = np.min(r_mag)
    r_max = np.max(r_mag)
    a = (r_max + r_min)/2
    # angle for period
    phi = np.arctan2(r_vec[:,1], r_vec[:,0])
    phi_un = np.unwrap(phi)
    delta_phi = phi_un[-1] - phi_un[0]
    cycles = delta_phi/(2*np.pi)
    T = (t[-1] - t[0])/cycles
    return T, a, r_min, r_max, r_mag

T, a, r_min, r_max, r_mag = compute_period_and_a(t, r1, r2)
# ensure T positive
T = abs(T)
# compute velocity via central difference for interior points
dt = t[2:] - t[:-2]
v1 = (r1[2:] - r1[:-2]) / dt[:,None]
v2 = (r2[2:] - r2[:-2]) / dt[:,None]
# relative positions at interior i
r1_mid = r1[1:-1]
r2_mid = r2[1:-1]
r_mid = r2_mid - r1_mid
r_mid_mag = np.linalg.norm(r_mid, axis=1)
# speeds
i1 = np.linalg.norm(v1, axis=1)
i2 = np.linalg.norm(v2, axis=1)
ratio_arr = i2 / i1
mass_ratio = np.median(ratio_arr)
# total mass
g = 6.67430e-11
M_total = 4*np.pi**2 * a**3 / (g * T**2)
# individual masses
m1 = M_total * mass_ratio/(1+mass_ratio)
m2 = M_total/(1+mass_ratio)
# compute instantaneous total energy at interior points
e_inst = 0.5*m1*(i1**2) + 0.5*m2*(i2**2) - g*m1*m2/r_mid_mag
# summary
print("T (s):", T)
print("a (m):", a)
print("M_total (kg):", M_total)
print("mass_ratio m1/m2:", mass_ratio)
print("m1 (kg):", m1)
print("m2 (kg):", m2)
print("Instantaneous E stats (J): min, max, median:", np.min(e_inst), np.max(e_inst), np.median(e_inst))
print("Analytic E (J):", -g*m1*m2/(2*a))

# Recompute instantaneous energy with COM subtraction
import numpy as np
# Data
t = df['time'].values
r1 = df[['star1_x','star1_y','star1_z']].values
t = df['time'].values
r2 = df[['star2_x','star2_y','star2_z']].values
# orbital elements from before
T, a, r_min, r_max, r_mag = compute_period_and_a(t, r1, r2)
T = abs(T)
g = 6.67430e-11
M_total = 4*np.pi**2 * a**3 / (g * T**2)
# compute v1,v2 at time indices 1 to n-2 (central difference)
dt = t[2:] - t[:-2]
v1 = (r1[2:] - r1[:-2]) / dt[:,None]
v2 = (r2[2:] - r2[:-2]) / dt[:,None]
# compute mass ratio from v-rel speeds
i1 = np.linalg.norm(v1, axis=1)
i2 = np.linalg.norm(v2, axis=1)
ratio_arr = i2 / i1
mass_ratio = np.median(ratio_arr)
# masses
q = mass_ratio
m1 = M_total * q/(1+q)
m2 = M_total/(1+q)
# COM velocity
v_cm = (m1 * v1 + m2 * v2) / (m1 + m2)
# relative velocities
v1_rel = v1 - v_cm
v2_rel = v2 - v_cm
# mid positions and r12 mid
r1_mid = r1[1:-1]
r2_mid = r2[1:-1]
r_mid = r2_mid - r1_mid
r_mid_mag = np.linalg.norm(r_mid, axis=1)
# KE and U
KE1 = 0.5 * m1 * np.sum(v1_rel**2, axis=1)
KE2 = 0.5 * m2 * np.sum(v2_rel**2, axis=1)
PE = -g * m1 * m2 / r_mid_mag
E_inst = KE1 + KE2 + PE
print('E_inst stats: min, max, mean, median')
print(np.min(E_inst), np.max(E_inst), np.mean(E_inst), np.median(E_inst))
print('Analytic E:', -g * m1 * m2 / (2*a))

# compute E using relative motion approach
import numpy as np
t = df['time'].values
r1 = df[['star1_x','star1_y','star1_z']].values
r2 = df[['star2_x','star2_y','star2_z']].values
# Derive masses and a, etc. as before
from numpy import unwrap, arctan2
# compute r and phi
r_vec = r2 - r1
r_mag = np.linalg.norm(r_vec, axis=1)
phi = np.arctan2(r_vec[:,1], r_vec[:,0])
phi_un = np.unwrap(phi)
delta_phi = phi_un[-1] - phi_un[0]
cycles = delta_phi/(2*np.pi)
T = abs((t[-1] - t[0])/cycles)
r_min = r_mag.min()
r_max = r_mag.max()
a = (r_min + r_max)/2
g = 6.67430e-11
M_total = 4*np.pi**2 * a**3 / (g * T**2)
# central diff v1, v2
dt = t[2:] - t[:-2]
v1 = (r1[2:] - r1[:-2]) / dt[:,None]
v2 = (r2[2:] - r2[:-2]) / dt[:,None]
# relative velocity
v_rel = v2 - v1
# r at interior
theta = r_vec[1:-1]
r_mid = r2[1:-1] - r1[1:-1]
r_mid_mag = np.linalg.norm(r_mid, axis=1)
# mass ratio q
i1 = np.linalg.norm(v1, axis=1)
i2 = np.linalg.norm(v2, axis=1)
ratio_arr = i2 / i1
q = np.median(ratio_arr)
m1 = M_total * q/(1+q)
m2 = M_total/(1+q)
mu = m1 * m2 / (m1 + m2)
# compute E_rel at interior points
E_rel = 0.5 * mu * np.sum(v_rel**2, axis=1) - g*m1*m2 / r_mid_mag
print('E_rel stats: min, max, mean, median')
print(np.min(E_rel), np.max(E_rel), np.mean(E_rel), np.median(E_rel))
print('Analytic E:', -g*m1*m2/(2*a))

{'answer': -4.062566794302988e+39}

import pandas as pd;
print(df.columns)
print(df.head())

# Compute time span
import numpy as np
time_span = df['time'].max() - df['time'].min()
print(time_span, len(df))

import numpy as np
# Extract data
t = df['time'].values
x1 = df['star1_x'].values
y1 = df['star1_y'].values
x2 = df['star2_x'].values
y2 = df['star2_y'].values
# relative position (vector from star2 to star1)
rx = x1 - x2
ry = y1 - y2
r_rel = np.sqrt(rx**2 + ry**2)
# Summary of r_rel
r_min = np.min(r_rel)
r_max = np.max(r_rel)
r_mean = np.mean(r_rel)
ecc = (r_max - r_min)/(r_max + r_min)
print("r_min=", r_min, "r_max=", r_max)
print("r_mean=", r_mean)
print("ecc estimate=", ecc)

# compute ratio r2/r1 for first and some other time points
idxs = [0, 1000, 2000, 3000, 4000, 5000, 6500-1]
ratios = []
r1 = np.sqrt(x1**2 + y1**2)
r2 = np.sqrt(x2**2 + y2**2)
for idx in idxs:
    ratios.append(r2[idx]/r1[idx])
print(ratios)

import numpy as np
import pandas as pd
# alias arrays
t = df['time'].values
r1 = np.vstack((df['star1_x'].values, df['star1_y'].values)).T  # shape (N,2)
r2 = np.vstack((df['star2_x'].values, df['star2_y'].values)).T
# relative positions
r_rel = r1 - r2  # shape (N,2)
r_rel_mag = np.linalg.norm(r_rel, axis=1)
# find pericenters: local minima of r_rel_mag
def find_local_minima(x):
    # return indices i where x[i] < x[i-1] and < x[i+1]
    # exclude first and last
    diff1 = x[1:-1] < x[:-2]
    diff2 = x[1:-1] < x[2:]
    return np.where(diff1 & diff2)[0] + 1
min_idxs = find_local_minima(r_rel_mag)
# get times of pericenters, ensure sorted
times_peri = t[min_idxs]
# print number of peri events and first few
overprint = (len(times_peri), times_peri[:10])
print("Number of pericenters (local minima):", overprint[0])
print("First few pericenter times:", overprint[1])
# compute period as differences between successive peri times
if len(times_peri) >= 2:
    period_diffs = np.diff(times_peri)
    P = np.mean(period_diffs)
else:
    P = t[-1] - t[0]
print("Estimated orbital period P (s):", P)
# compute r_min, r_max, a_rel, e
r_min = np.min(r_rel_mag)
r_max = np.max(r_rel_mag)
a_rel = (r_min + r_max) / 2
ecc = (r_max - r_min) / (r_max + r_min)
print("r_min (m):", r_min)
print("r_max (m):", r_max)
print("a_rel (semi-major axis) (m):", a_rel)
print("Eccentricity estimate:", ecc)
# gravitational constant
G = 6.67430e-11
# compute total mass M_total from Kepler's third law
M_total = 4 * np.pi**2 * a_rel**3 / (G * P**2)
print("Estimated total mass M (kg):", M_total)

# detect local maxima
x = r_rel_mag
diff1 = x[1:-1] > x[:-2]
diff2 = x[1:-1] > x[2:]
max_idxs = np.where(diff1 & diff2)[0] + 1
print(len(max_idxs), max_idxs[:10])
# detect endpoints maxima
end_max = []
if x[0] > x[1]: end_max.append(0)
if x[-1] > x[-2]: end_max.append(len(x)-1)
print("end_max", end_max)

# Compute sum x1+x2, y1+y2
sx = x1 + x2
y_sum = y1 + y2
print("x1+x2 min, max:", np.min(sx), np.max(sx))
print("y1+y2 min, max:", np.min(y_sum), np.max(y_sum))
# Also compute differences dx= x1* M1 + x2*M2 for M1=M2=1
dx0 = sx
# Variation range
print("x1+x2 variation:", np.max(sx)-np.min(sx))
print("y1+y2 variation:", np.max(y_sum)-np.min(y_sum))

print(r_rel_mag[0], r_rel_mag[-1])

We have a clean 2‐D data set of two stars moving in the plane (z≡0).  Our ultimate goal is the total mechanical energy 

    E = K + U  
      = ½ M₁ |v₁|² + ½ M₂ |v₂|²  –  G M₁ M₂ / |r₁–r₂|  

But as everyone knows, for an isolated two–body system one can equally well work in the relative coordinate r = r₁–r₂ with the reduced mass μ = M₁ M₂/(M₁+M₂).  The conserved energy becomes

    E = ½ μ |v_rel|²  –  G M₁ M₂ / |r_rel|

and, for a bound ellipse, one also has the closed‐form

    E  =  –  G M₁ M₂  
             ───────────  
           2 a_rel  

where a_rel is the semi‐major axis of the relative orbit.

So our recipe is:

1.  From the raw (x₁,y₁), (x₂,y₂) & t data  
    a) build the relative separation r_rel(t) = r₁(t) – r₂(t) and measure its pericenter and apocenter to get  
       • r_min, r_max ⇒ a_rel = (r_min + r_max)/2  
       • eccentricity e = (r_max – r_min)/(r_max + r_min)  
       • orbital period P from successive pericenter (or apocenter) timings.  
    b) Use Kepler’s 3rd law in SI:  
       
         M_total ≡ M₁ + M₂  =  4 π² a_rel³  
                              ───────────  
                                 G P²  

2.  To split M_total into M₁ and M₂ we need the mass ratio q = M₂/M₁.  
    •  In the barycentric frame r₁_center = +[M₂/(M₁+M₂)] r_rel and  
                        r₂_center = –[M₁/(M₁+M₂)] r_rel  
       ⇒  |r₂_center|/|r₁_center| = M₁/M₂ = 1/q is *constant*.  
    •  Our data are not originally in the barycenter frame, so we must (numerically) solve for q by demanding that the inferred center‐of‐mass  

           R_cm(t; q) = [M₁ r₁(t) + M₂ r₂(t)]/(M₁+M₂)

       moves on a straight line (i.e. constant velocity) with minimal residual scatter.  
    •  Once q is found, we have:  
          
           M₁ =  M_total  / (1 + q)  
           M₂ =  q M₁  

3.  With M₁, M₂ and a_rel in hand, plug into  

       E = –  G M₁ M₂  
               ───────────  
             2 a_rel  

   to get the conserved total energy in Joules.  

Alternatively (and as a consistency check) we can transform into the barycentric frame (subtract off R_cm and its velocity), finite–difference the positions to get barycentric velocities v₁, v₂, and then at any epoch form  

    K = ½ M₁ |v₁|² + ½ M₂ |v₂|²  
    U = –G M₁ M₂ / |r_rel|  
    E_check = K + U  

and verify it sits on the same constant value we get from the –G M₁ M₂/(2a) formula.  

Next steps:  
– Identify pericenter/apocenter to get P, r_min, r_max, a_rel  
– Compute M_total = 4π²a_rel³/(G P²)  
– Numerically solve for q by demanding R_cm(t;q) is best‐fit by a straight line  
– Compute M₁, M₂, and E_total = –G M₁ M₂/(2 a_rel)  

I will now move to the code‐interpreter and implement each block in turn.

import numpy as np
# 1. Compute r_rel and r_min, r_max
t = df['time'].values
r1 = np.vstack((df['star1_x'], df['star1_y'])).T
r2 = np.vstack((df['star2_x'], df['star2_y'])).T
r_rel = r1 - r2
r_mag = np.linalg.norm(r_rel, axis=1)
r_min = r_mag.min()
r_max = r_mag.max()
a_rel = 0.5*(r_min + r_max)
ecc = (r_max - r_min)/(r_max + r_min)
print(f"r_min = {r_min:.3e} m")
print(f"r_max = {r_max:.3e} m")
print(f"a_rel = {a_rel:.3e} m")
print(f"e = {ecc:.6f}")

# 2. Compute P and M_total
G = 6.67430e-11
t0, tN = t[0], t[-1]
time_span = tN - t0
P = 2 * time_span  # assume from pericenter at t0 to apocenter at tN
M_total = 4 * np.pi**2 * a_rel**3 / (G * P**2)
print(f"t0 = {t0:.3e} s, tN = {tN:.3e} s, time_span = {time_span:.3e} s")
print(f"Estimated orbital period P = {P:.3e} s = {P/ (3600*24*365):.3f} years")
print(f"Estimated total mass M_total = {M_total:.3e} kg = {M_total/1.98847e30:.3f} M_sun") # solar mass ~1.98847e30 kg

import numpy as np
import pandas as pd
# Prepare data arrays
t = df['time'].values
tau = t - t[0]  # time from first measurement
x1 = df['star1_x'].values
y1 = df['star1_y'].values
x2 = df['star2_x'].values
y2 = df['star2_y'].values

# SSR(q) function: residual sum of squares of COM linear fit

def SSR(q):
    M_ratio = q
    # COM positions
    com_x = (x1 + M_ratio * x2) / (1 + M_ratio)
    com_y = (y1 + M_ratio * y2) / (1 + M_ratio)
    # Fit com_x ~ a_x + b_x * tau
    # Solve for b_x = cov(com_x, tau)/var(tau)
    mean_tau = np.mean(tau)
    mean_com_x = np.mean(com_x)
    mean_com_y = np.mean(com_y)
    # cov(x,tau)
    cov_x_tau = np.mean((com_x - mean_com_x) * (tau - mean_tau))
    var_tau = np.mean((tau - mean_tau)**2)
    b_x = cov_x_tau / var_tau
    a_x = mean_com_x - b_x * mean_tau
    # same for y
    cov_y_tau = np.mean((com_y - mean_com_y) * (tau - mean_tau))
    b_y = cov_y_tau / var_tau
    a_y = mean_com_y - b_y * mean_tau
    # residuals
    rx = com_x - (a_x + b_x * tau)
    ry = com_y - (a_y + b_y * tau)
    return np.sum(rx*rx + ry*ry)

# Quick test SSR at some q values
qs = [0.1, 0.5, 1, 2, 5, 10]
ssr_vals = [(q, SSR(q)) for q in qs]
print(ssr_vals)

# find q minimizing SSR with brute log sampling
grid = np.logspace(-1, 1, 200)  # from 0.1 to 10
ssr_grid = [SSR(q) for q in grid]
min_idx = np.argmin(ssr_grid)
q_best = grid[min_idx]
ssr_best = ssr_grid[min_idx]
print("Best q in [0.1,10]:", q_best, "SSR=", ssr_best)
# refine around q_best
grid2 = np.linspace(max(q_best*0.5, 1e-4), q_best*1.5, 200)
ssr2 = [SSR(q) for q in grid2]
min2 = np.argmin(ssr2)
q_best2 = grid2[min2]
print("Refined q in [{:.3e},{:.3e}]: bruted says {:.6f}".format(grid2[0],grid2[-1],q_best2))

from scipy.optimize import minimize_scalar
# find optimal q via bounded minimization
do_bounds = (0.2, 2.0)
res = minimize_scalar(SSR, bounds=do_bounds, method='bounded', options={'xatol':1e-6})
q_opt = res.x
print(res)

# compute M1, M2 and energy via formula
q = q_opt  # from SSR minimization
M_total = M_total  # from above
M1 = M_total / (1 + q)
M2 = M_total - M1
# total energy via E = -G M1 M2 / (2 a_rel)
E_formula = - G * M1 * M2 / (2 * a_rel)
print(f"q = {q:.6f}")
print(f"M1 = {M1:.3e} kg = {M1/1.98847e30:.3f} M_sun")
print(f"M2 = {M2:.3e} kg = {M2/1.98847e30:.3f} M_sun")
print(f"Total energy (via -G M1 M2/(2a)): E = {E_formula:.3e} J")
# Now approximate E(t) = K + U across sample times
dt = np.diff(t)
# velocities via central differences
v1 = np.zeros_like(r1)
v2 = np.zeros_like(r2)
# central for interior
i0, iN = 1, len(t)-1
for i in range(1, len(t)-1):
    dt_c = t[i+1] - t[i-1]
    v1[i] = (r1[i+1] - r1[i-1]) / dt_c
    v2[i] = (r2[i+1] - r2[i-1]) / dt_c
# endpoints forward/backward
dt_f = t[1] - t[0]
v1[0] = (r1[1] - r1[0]) / dt_f
v2[0] = (r2[1] - r2[0]) / dt_f
dt_b = t[-1] - t[-2]
v1[-1] = (r1[-1] - r1[-2]) / dt_b
v2[-1] = (r2[-1] - r2[-2]) / dt_b
# K(t)
K = 0.5 * M1 * np.sum(v1**2, axis=1) + 0.5 * M2 * np.sum(v2**2, axis=1)
# U(t)
U = - G * M1 * M2 / r_mag
E_t = K + U
print(f"E(t) stats: min {E_t.min():.3e}, max {E_t.max():.3e}, mean {np.mean(E_t):.3e}")

import numpy as np
from scipy.optimize import minimize_scalar
# load arrays
t = df['time'].values
r1 = np.vstack((df['star1_x'].values, df['star1_y'].values)).T
r2 = np.vstack((df['star2_x'].values, df['star2_y'].values)).T
# 1) r_min, r_max, a
r_rel = r1 - r2
r_mag = np.linalg.norm(r_rel, axis=1)
r_min = r_mag.min()
r_max = r_mag.max()
a_rel = 0.5*(r_min + r_max)
# 2) period P = 2*(t[-1] - t[0])
t0 = t[0]	N = t[-1]
time_span = N - t0
P = 2 * time_span
# constants
g = 6.67430e-11
pi = np.pi
# total mass
epsilon = 4*pi**2 * a_rel**3 / (g * P**2)
# SSR function for q
x1, y1 = r1[:,0], r1[:,1]
x2, y2 = r2[:,0], r2[:,1]
tau = t - t0
mean_tau = np.mean(tau)
var_tau = np.mean((tau - mean_tau)**2)

def SSR(q):
    com_x = (x1 + q*x2)/(1+q)
    com_y = (y1 + q*y2)/(1+q)
    # linear fit coefficients
    cov_x = np.mean((com_x - com_x.mean()) * (tau - mean_tau))
    cov_y = np.mean((com_y - com_y.mean()) * (tau - mean_tau))
    bx = cov_x/var_tau
    by = cov_y/var_tau
    ax = com_x.mean() - bx*mean_tau
    ay = com_y.mean() - by*mean_tau
    rx = com_x - (ax + bx*tau)
    ry = com_y - (ay + by*tau)
    return np.sum(rx*rx + ry*ry)
# minimize SSR
es = minimize_scalar(SSR, bounds=(0.1,10), method='bounded', options={'xatol':1e-6})
q_opt = res.x
# Masses
M_total = epsilon
epsilon_Mtot = M_total
M1 = M_total/(1+q_opt)
M2 = M_total - M1
# energy
e_formula = - g * M1 * M2 / (2 * a_rel)
# print results
print(r_min, r_max, a_rel)
print(P)
print("M_total=", M_total)
print("q_opt=", q_opt)
print("M1=", M1)
print("M2=", M2)
print("E=", e_formula)

import numpy as np
# reload data from df
t = df['time'].values
t0 = t[0]
# positions
r1 = np.vstack((df['star1_x'], df['star1_y'])).T
r2 = np.vstack((df['star2_x'], df['star2_y'])).T
# relative separation vector and magnitude
r_rel = r1 - r2
r_mag = np.linalg.norm(r_rel, axis=1)
# find pericenter and apocenter indices
i_min = np.argmin(r_mag)
i_max = np.argmax(r_mag)
# times at pericenter and apocenter
t_min = t[i_min]
t_max = t[i_max]
# period estimate
time_half = abs(t_max - t_min)
P = 2 * time_half
# semi-major axis
a_rel = 0.5 * (r_mag[i_min] + r_mag[i_max])
# G constant
G = 6.67430e-11
pi = np.pi
# total mass
M_total = 4 * pi**2 * a_rel**3 / (G * P**2)
# print details
def sci(x): return f"{x:.6e}"
print("i_min (peri), i_max (apo):", i_min, i_max)
print("r_min (m)    =", sci(r_mag[i_min]))
print("r_max (m)    =", sci(r_mag[i_max]))
print("a_rel (m)    =", sci(a_rel))
print("t_min (s)    =", sci(t_min))
print("t_max (s)    =", sci(t_max))
print("time half    =", sci(time_half), "s")
print("P estimate   =", sci(P), "s =", P/(3600*24*365), "yr")
print("M_total      =", sci(M_total), "kg =", M_total/1.98847e30, "M_sol")

{'answer': -2.184e+42}

# Inspect the DataFrame
def print_df_info():
    print("Columns:", df.columns.tolist())
    print(df.head())

print_df_info()

# Check loaded global variables (non-private)
print([name for name in globals() if not name.startswith('_')])

print(len(df))

# Compute relative positions
r = df[['star2_x', 'star2_y', 'star2_z']].values - df[['star1_x', 'star1_y', 'star1_z']].values
mean_rx = r[:,0].mean()
mean_ry = r[:,1].mean()
mean_rz = r[:,2].mean()
print(mean_rx, mean_ry, mean_rz)

import numpy as np
# Extract positions & times
t = df['time'].values
r1 = df[['star1_x','star1_y','star1_z']].values
r2 = df[['star2_x','star2_y','star2_z']].values
# Time step differences
dt1 = np.diff(t)
print('dt unique values (first 10):', np.unique(np.round(dt1, 10))[:10])
print('All dt equal? ', np.allclose(dt1, dt1[0]))
# Use central difference for acceleration
dt = dt1[0]
n = len(df)
a1 = (r1[2:] - 2*r1[1:-1] + r1[:-2])/(dt**2)
a2 = (r2[2:] - 2*r2[1:-1] + r2[:-2])/(dt**2)
r = r2[1:-1] - r1[1:-1]
r_mag = np.linalg.norm(r, axis=1)
a1_mag = np.linalg.norm(a1, axis=1)
a2_mag = np.linalg.norm(a2, axis=1)
# gravitational constant
G = 6.67430e-11
# Compute mass estimates
tiny = 1e-20
m2_est = a1_mag * (r_mag**2) / G
m1_est = a2_mag * (r_mag**2) / G
# Check nan or inf
print('m2_est stats:', np.nanmin(m2_est), np.nanmedian(m2_est), np.nanmax(m2_est))
print('m1_est stats:', np.nanmin(m1_est), np.nanmedian(m1_est), np.nanmax(m1_est))

import numpy as np
G = 6.67430e-11
# Data
import numpy as np
# loaded: df
t = df['time'].values
r1 = df[['star1_x','star1_y','star1_z']].values
r2 = df[['star2_x','star2_y','star2_z']].values
# dt (assume constant)
dt = np.mean(np.diff(t))
# acceleration central diff
a1 = (r1[2:] - 2*r1[1:-1] + r1[:-2]) / (dt**2)
a2 = (r2[2:] - 2*r2[1:-1] + r2[:-2]) / (dt**2)
# r vector at 1:-1
r = r2[1:-1] - r1[1:-1]
r_mag = np.linalg.norm(r, axis=1)
# masses
m2_est = np.linalg.norm(a1,axis=1)* (r_mag**2)/G
m1_est = np.linalg.norm(a2,axis=1)* (r_mag**2)/G
m1 = np.median(m1_est)
m2 = np.median(m2_est)
print("Estimated masses (median): m1=",m1,"kg, m2=",m2,"kg")
# velocities central diff
v1 = (r1[2:] - r1[:-2])/(2*dt)
v2 = (r2[2:] - r2[:-2])/(2*dt)
# kinetic energies
K1 = 0.5*m1 * np.sum(v1**2,axis=1)
K2 = 0.5*m2 * np.sum(v2**2,axis=1)
K = K1+K2
# potential energies
U = -G*m1*m2/r_mag
# total energy
E = K+U
print("K stats: min, mean, max =", np.min(K), np.mean(K), np.max(K))
print("U stats: min, mean, max =", np.min(U), np.mean(U), np.max(U))
print("E stats: min, mean, max, std =", np.min(E), np.mean(E), np.max(E), np.std(E))
print("Mean E:", np.mean(E))
print("Median E:", np.median(E))

# compute r_mag
d = df[['star2_x','star2_y','star2_z']].values - df[['star1_x','star1_y','star1_z']].values
rmag = np.linalg.norm(d,axis=1)
print('r_mag stats: min, mean, max, std:', np.min(rmag), np.mean(rmag), np.max(rmag), np.std(rmag))
# find extremal indices
imin = np.argmin(rmag)
imax = np.argmax(rmag)
print('min index, time, rmag:', imin, df['time'].iloc[imin], rmag[imin])
print('max index, time, rmag:', imax, df['time'].iloc[imax], rmag[imax])

import numpy as np
# Data
import numpy as np
t = df['time'].values
r1 = df[['star1_x','star1_y','star1_z']].values
r2 = df[['star2_x','star2_y','star2_z']].values
# time step
dt = np.mean(np.diff(t))
# velocities
v1 = (r1[2:] - r1[:-2])/(2*dt)
v2 = (r2[2:] - r2[:-2])/(2*dt)
# accelerations
a1 = (r1[2:] - 2*r1[1:-1] + r1[:-2])/(dt**2)
a2 = (r2[2:] - 2*r2[1:-1] + r2[:-2])/(dt**2)
# relative positions
r = r2[1:-1] - r1[1:-1]
r_mag = np.linalg.norm(r,axis=1)
# gravitational constant
G = 6.67430e-11
# mass estimates
m2_est = np.linalg.norm(a1,axis=1) * (r_mag**2) / G
m1_est = np.linalg.norm(a2,axis=1) * (r_mag**2) / G
m1 = np.median(m1_est)
m2 = np.median(m2_est)
# kinetic energy
e1 = 0.5*m1*np.sum(v1**2,axis=1)
e2 = 0.5*m2*np.sum(v2**2,axis=1)
K = e1 + e2
# potential energy
U = -G*m1*m2 / r_mag
e = K + U
# statistics
e_mean = np.mean(e)
e_min = np.min(e)
e_max = np.max(e)
e_std = np.std(e)
print('m1, m2 =', m1, m2)
print('E mean =', e_mean)
print('E median =', np.median(e))
print('E min, max, std =', e_min, e_max, e_std)

{'answer': 1.3771392221670886e+40}

import numpy as np
import pandas as pd
# Inspect the DataFrame
print(df.head())
print('\nDataFrame info:')
print(df.info())
print('\nSummary statistics:')
print(df.describe())

print(df[['star1_x','star1_y','star2_x','star2_y']].head())

import numpy as np
# Extract arrays
t = df['time'].values
r1 = df[['star1_x','star1_y']].values
r2 = df[['star2_x','star2_y']].values
# Relative vector
r_rel = r2 - r1
# Magnitude of relative vector
r_mag = np.sqrt((r_rel[:,0])**2 + (r_rel[:,1])**2)
# Angle and unwrap for period
theta = np.arctan2(r_rel[:,1], r_rel[:,0])
unwrap_theta = np.unwrap(theta)
# Slope of angle vs time
slope = (unwrap_theta[-1] - unwrap_theta[0])/(t[-1] - t[0])
# Period
T = 2*np.pi/slope
# Radial extents and semimajor axis estimation
r_min = np.min(r_mag)
r_max = np.max(r_mag)
a = 0.5*(r_min + r_max)
# Compute total mass from Kepler's 3rd law
G = 6.674e-11
M_tot = 4*np.pi**2 * a**3/(G * T**2)
# Compute mass ratio parameter alpha = m1/M_tot by minimising COM variance
# COM_i = alpha*r1_i + (1-alpha)*r2_i = r2_i + alpha*(r1_i - r2_i)
d = r1 - r2
dx = d[:,0]
dy = d[:,1]
r2x = r2[:,0]
r2y = r2[:,1]
# Means
Edx = np.mean(dx)
Edy = np.mean(dy)
Er2x = np.mean(r2x)
Er2y = np.mean(r2y)
# Variances of d
var_dx = np.mean(dx*dx) - Edx**2
var_dy = np.mean(dy*dy) - Edy**2
# Covariances
Cov_x = np.mean(r2x*dx) - Er2x*Edx
Cov_y = np.mean(r2y*dy) - Er2y*Edy
alpha = -(Cov_x + Cov_y)/(var_dx + var_dy)
# Masses
m1 = alpha * M_tot
m2 = (1 - alpha) * M_tot
# Total energy (internal)
E_int = -G * m1 * m2/(2 * a)
# Print results
print("Period T =", T, "s")
print("r_min =", r_min, "m, r_max =", r_max, "m, a =", a, "m")
print("Total mass M_tot =", M_tot, "kg")
print("Mass fraction alpha = m1/M_tot =", alpha)
print("m1 =", m1, "kg, m2 =", m2, "kg, m1/m2 =", m1/m2)
print("Total energy E_int =", E_int, "J")

# Recalculate with corrections
import numpy as np
# Coordinates
t = df['time'].values
r1 = df[['star1_x','star1_y','star1_z']].values
r2 = df[['star2_x','star2_y','star2_z']].values
# Relative vector and magnitude
r_rel = r2 - r1
r_mag = np.linalg.norm(r_rel, axis=1)
# Get theta and unwrap
theta = np.arctan2(r_rel[:,1], r_rel[:,0])
unwrap_theta = np.unwrap(theta)
# Use linear regression for unwrap_theta vs t
a = np.vstack([t, np.ones_like(t)]).T
slope, intercept = np.linalg.lstsq(a, unwrap_theta, rcond=None)[0]
# Period positive
tau = 2*np.pi/abs(slope)
# Semimajor axis
r_min = np.min(r_mag)
r_max = np.max(r_mag)
semimajor = 0.5*(r_min+r_max)
# total mass from Kepler's third
G = 6.674e-11
M_tot = 4*np.pi**2 * semimajor**3/(G * tau**2)
# Determine alpha by minimizing variance of COM
# d = r1 - r2
d = r1[:, :2] - r2[:, :2]  # use x,y
r2_xy = r2[:, :2]
# Compute covariance and variance components
Ed = d.mean(axis=0)
Er2 = r2_xy.mean(axis=0)
var_d = np.mean((d - Ed)**2)
Cov_r2_d = np.mean((r2_xy - Er2) * (d - Ed))  # elementwise multiply then mean sum across axis?
# Cov_r2_d should be sum of dot products? But we need Cov(r2x,dx)+Cov(r2y,dy)
Cov_x = np.mean((r2_xy[:,0]-Er2[0]) * (d[:,0]-Ed[0]))
Cov_y = np.mean((r2_xy[:,1]-Er2[1]) * (d[:,1]-Ed[1]))
Cov_sum = Cov_x + Cov_y
var_sum = ((d[:,0]-Ed[0])**2 + (d[:,1]-Ed[1])**2).mean()
alpha = -Cov_sum/var_sum
# masses
global_m1 = alpha * M_tot
global_m2 = (1-alpha) * M_tot
# total energy
E_total = -G*global_m1*global_m2/(2*semimajor)
# Print
debug = {
    'slope': slope,
    'tau': tau,
    'a': semimajor,
    'M_tot': M_tot,
    'alpha': alpha,
    'm1': global_m1,
    'm2': global_m2,
    'energy': E_total,
}
for k,v in debug.items():
    print(f"{k} = {v}")

import numpy as np
# use all existing df
# Coordinates
t = df['time'].values
r1 = df[['star1_x','star1_y','star1_z']].values
r2 = df[['star2_x','star2_y','star2_z']].values
# Relative vector and magnitude
r_rel = r2 - r1
r_mag = np.linalg.norm(r_rel, axis=1)
# Get semimajor axis
r_min = r_mag.min()
 r_max = r_mag.max()

import numpy as np
import pandas as pd
# Use df loaded
t = df['time'].values
r1 = df[['star1_x','star1_y']].values
r2 = df[['star2_x','star2_y']].values
# relative vector and magnitude
r_rel = r2 - r1
r_mag = np.linalg.norm(r_rel, axis=1)
# semimajor axis
a = (r_mag.max() + r_mag.min()) / 2.0
# orbital angle
theta = np.unwrap(np.arctan2(r_rel[:,1], r_rel[:,0]))
# linear regression for slope
a_coef, b_coef = np.polyfit(t, theta, 1)
slope = a_coef
T = 2 * np.pi / abs(slope)
# total mass M_tot
g = 6.67408e-11
M_tot = 4 * np.pi**2 * a**3 / (g * T**2)
# get alpha via COM variance minimization
d = r1 - r2
# covariances and variances
dx = d[:,0];
dy = d[:,1]
# r2 coords
y2 = r2[:,1];
x2 = r2[:,0]
# Cov_x = Cov(x2,dx)
Cov_x = np.mean((x2 - x2.mean()) * (dx - dx.mean()))
Cov_y = np.mean((y2 - y2.mean()) * (dy - dy.mean()))
cov_sum = Cov_x + Cov_y
# variances d
var_dx = np.mean((dx - dx.mean())**2)
var_dy = np.mean((dy - dy.mean())**2)
var_sum = var_dx + var_dy
alpha = - cov_sum / var_sum
# masses
m1 = alpha * M_tot
m2 = M_tot - m1
# energy
E_tot = - g * m1 * m2 / (2 * a)

# print results
print("a =", a)
print("T =", T)
print("M_tot =", M_tot)
print("alpha =", alpha)
print("m1 =", m1)
print("m2 =", m2)
print("Energy =", E_tot)

# Cross-check energy from KE+U at one snapshot
import numpy as np
# using m1, m2, G, df
g = 6.67408e-11
# parameters from above
a = 3494100967130.564
T = 723599523.003957
M_tot = 4.819235352080961e+31
alpha = 0.8729508196725002
m1 = alpha*M_tot
m2 = M_tot - m1
# choose index near middle
i = len(df)//2
# times t[i-1], t[i], t[i+1]
t = df['time'].values
# positions
pos1 = df[['star1_x','star1_y','star1_z']].values
pos2 = df[['star2_x','star2_y','star2_z']].values
# approximate velocities at i by central diff
dt_forward = t[i+1] - t[i]
dt_backward = t[i] - t[i-1]
# central velocity average
v1_forward = (pos1[i+1] - pos1[i])/dt_forward
v1_backward = (pos1[i] - pos1[i-1])/dt_backward
v1 = 0.5*(v1_forward + v1_backward)

v2_forward = (pos2[i+1] - pos2[i])/dt_forward
v2_backward = (pos2[i] - pos2[i-1])/dt_backward
v2 = 0.5*(v2_forward + v2_backward)
# COM velocity v_COM = (m1*v1 + m2*v2)/(m1+m2)
v_COM = (m1*v1 + m2*v2)/M_tot
# relative velocity with COM frame v_rel1 = v1 - v_COM, v_rel2 = v2 - v_COM, but relative velocity between masses? Could use v_rel = v2 - v1
# Actually relative velocity is difference in velocities in inertial: v_rel = v2 - v1
v_rel = v2 - v1
# reduced mass mu = m1*m2/M_tot
mu = m1*m2/M_tot
# KE internal = 0.5*mu*|v_rel|^2
KE_int = 0.5*mu*np.dot(v_rel,v_rel)
# potential U = -G*m1*m2/|r2-r1|
r_rel = pos2[i] - pos1[i]
U = -g*m1*m2/np.linalg.norm(r_rel)
# internal energy E = KE+U
E_int_snapshot = KE_int + U
print("m1, m2", m1, m2)
print("v1", v1)
print("v2", v2)
print("v_rel", v_rel)
print("KE_int", KE_int)
print("U", U)
print("E_int_snapshot", E_int_snapshot)
print("E_expected", -g*m1*m2/(2*a))

{'answer': -2.46004968573789e+39}

import pandas as pd
# Inspect columns and head of the DataFrame
df.head()
print(df.columns)
print(df.info())

# Display some sample rows of positions to see values magnitude and pattern
print(df[['star1_x','star1_y','star1_z','star2_x','star2_y','star2_z']].head(10))

# Compute mean positions
mean_star1 = df[['star1_x','star1_y','star1_z']].mean()
mean_star2 = df[['star2_x','star2_y','star2_z']].mean()
print(mean_star1, mean_star2)

# Check time differences
dt = df['time'].iloc[1] - df['time'].iloc[0]
print(dt, df['time'].min(), df['time'].max())

# Compute COM position under equal masses assumption as average of star1 and star2 positions
# Actually COM_x = (star1_x + star2_x)/2
df['com_x_eq'] = (df['star1_x'] + df['star2_x']) / 2
df['com_y_eq'] = (df['star1_y'] + df['star2_y']) / 2
df['com_z_eq'] = (df['star1_z'] + df['star2_z']) / 2
print(df[['com_x_eq','com_y_eq','com_z_eq']].head())
# Compute mean and std of COM positions
df[['com_x_eq','com_y_eq','com_z_eq']].describe()

import numpy as np
# Compute separation vector and its magnitude for each time
sep_vec = df[['star1_x','star1_y','star1_z']].values - df[['star2_x','star2_y','star2_z']].values
d_sep = np.linalg.norm(sep_vec, axis=1)
print('Separation mean:', d_sep.mean(), 'std:', d_sep.std(), 'min:', d_sep.min(), 'max:', d_sep.max())

print([var for var in globals().keys() if not var.startswith('_')])

# compute first 10 separations and print
dprint = d_sep[:20]
for i, val in enumerate(dprint): print(i, val)

df[['com_x_eq','com_y_eq','com_z_eq']].describe()

import numpy as np, pandas as pd
r = d_sep
t = df['time'].values
# find local minima indices
inds_min = np.where((r[1:-1] < r[:-2]) & (r[1:-1] < r[2:]))[0] + 1
# Also threshold: ensure local minima where r less than some percentile? But no, it's stable.
print('Local minima count:', len(inds_min))
# times and separation values at minima
t_min = t[inds_min]
r_min_vals = r[inds_min]
print('r_min_vals:', r_min_vals[:10])
print('First 5 r_min_vals and t_min:')
for i in range(min(5, len(inds_min))): print(i, r_min_vals[i], t_min[i])
# differences between successive minima
dt_min = np.diff(t_min)
print('dt_min:', dt_min[:10])
# statistics of dt_min
dt_stats = {'mean': dt_min.mean(), 'std': dt_min.std(), 'min': dt_min.min(), 'max': dt_min.max()}
print('dt_min stats:', dt_stats)
# local maxima indices
inds_max = np.where((r[1:-1] > r[:-2]) & (r[1:-1] > r[2:]))[0] + 1
print('Local maxima count:', len(inds_max))
t_max = t[inds_max]
r_max_vals = r[inds_max]
print('First 5 r_max_vals:', r_max_vals[:5])
print('t_max first 5 and r_max:')
for i in range(min(5, len(inds_max))): print(i, r_max_vals[i], t_max[i])

dt_max = np.diff(t_max)
print('dt_max stats:', {'mean': dt_max.mean(), 'std': dt_max.std(), 'min': dt_max.min(), 'max': dt_max.max()})
# compute average r_peri and r_apo
r_peri = np.mean(r_min_vals)
r_apo = np.mean(r_max_vals)
# compute a and e
a = (r_peri + r_apo) / 2
e = (r_apo - r_peri) / (r_apo + r_peri)
print('r_peri:', r_peri, 'r_apo:', r_apo, 'a:', a, 'e:', e)

# Print first and last COM positions assuming equal masses
print(df[['com_x_eq','com_y_eq','com_z_eq']].iloc[0])
print(df[['com_x_eq','com_y_eq','com_z_eq']].iloc[-1])

import numpy as np
# Load required arrays
t = df['time'].values
r1 = df[['star1_x','star1_y','star1_z']].values
r2 = df[['star2_x','star2_y','star2_z']].values
# Define objective: sum of variances of COM_x and COM_y

def com_variance(k):
    # COM coords
    com = (k * r1 + r2) / (k + 1)
    # We require COM stationary: minimize variance of x and y
    var_x = np.var(com[:,0])
    var_y = np.var(com[:,1])
    # var_z is zero
    return var_x + var_y
# Let's try some k values around 1
test_ks = np.logspace(-1,1,50)  # from 0.1 to 10
vals = [com_variance(k) for k in test_ks]
# find best k among test_ks
best_idx = np.argmin(vals)
best_k_rough = test_ks[best_idx]
print('Best k rough:', best_k_rough, 'variance:', vals[best_idx])
# refine around best_k_rough using golden section search
from scipy.optimize import minimize_scalar
res = minimize_scalar(com_variance, bounds=(test_ks[max(best_idx-1,0)], test_ks[min(best_idx+1,len(test_ks)-1)]),
                      method='bounded', options={'xatol':1e-6})
best_k = res.x
print('Optimization result:', res)
print('Best k refined:', best_k, 'variance:', com_variance(best_k))
# Check variance at k=1
print('Variance at k=1:', com_variance(1.0))
# Also at best_k: print range of COM positions
com_best = (best_k * r1 + r2) / (best_k + 1)
print('COM_x range:', com_best[:,0].max() - com_best[:,0].min(),
      'COM_y range:', com_best[:,1].max() - com_best[:,1].min())

# Compute masses and total energy
g = 6.67430e-11
# Use values from above computations
r_peri = 1.3975715184349148e+12  # from r_peri
r_apo  = 5.590630644020784e+12   # from r_apo
a = (r_peri + r_apo) / 2
P = 721067172.379743  # mean orbital period
k = 6.870967699549907  # mass ratio m1/m2
# Compute total mass
M_tot = 4 * np.pi**2 * a**3 / (g * P**2)
# Compute individual masses
m2 = M_tot / (k + 1)
m1 = M_tot - m2
# Compute total energy E = -G m1 m2 / (2a)
E = -g * m1 * m2 / (2 * a)
# Print results
print('a:', a)
print('P:', P)
print('M_tot:', M_tot, 'kg')
print('m1:', m1, 'kg')
print('m2:', m2, 'kg')
print('E_total:', E, 'J')

{'answer': -2.4947086166657293e+39}

import numpy as np
import pandas as pd
# Let's inspect the dataframe
df.head()

print(df.columns)
print(df.head())

print(df[['star2_x', 'star2_y', 'star2_z']].head())

print(len(df))
print(df['time'].min(), df['time'].max())
print(df['time'].max() - df['time'].min())

# Compute relative separation distances
r_vec = df[['star2_x', 'star2_y', 'star2_z']].values - df[['star1_x', 'star1_y', 'star1_z']].values
r = np.linalg.norm(r_vec, axis=1)
r_min = r.min()
r_max = r.max()
r_mean = r.mean()
r_min, r_max, r_max/r_min, r_mean, (r_max-r_min)/r_mean

print(r_min, r_max, r_max/r_min, r_mean, (r_max-r_min)/r_mean)

import numpy as np
# Compute r distances and times
times = df['time'].values
pos1 = df[['star1_x','star1_y','star1_z']].values
pos2 = df[['star2_x','star2_y','star2_z']].values
r = np.linalg.norm(pos2-pos1,axis=1)
# find local minima
def find_local_minima(r):
    minima = []
    for i in range(1,len(r)-1):
        if r[i] < r[i-1] and r[i] < r[i+1]:
            minima.append(i)
    return minima
minima_idx = find_local_minima(r)
# print number and some values
mins_r = [(i,r[i]) for i in minima_idx[:10]]
print("Total local minima detected:", len(minima_idx))
print("First 10 minima indices and r:", mins_r)
# times at minima
times_minima = times[minima_idx]
# compute diffs
diffs = np.diff(times_minima)
print("First 10 diffs:", diffs[:10])
print("All diffs shape:", diffs.shape)
# filter diffs >1e9 s (approx >30 years) 
deltas = diffs[diffs>1e9]
print("Filtered diffs >1e9 s, count:", len(deltas))
print("Filtered diffs sample:", deltas[:10])
print("Mean period (s):", np.mean(deltas))
print("Std of periods (s):", np.std(deltas))

import numpy as np
# positions from df
times = df['time'].values
pos1 = df[['star1_x','star1_y','star1_z']].values
pos2 = df[['star2_x','star2_y','star2_z']].values
r = np.linalg.norm(pos2-pos1,axis=1)
# find local minima (periapsis)
def find_local_minima(r):
    minima = []
    for i in range(1,len(r)-1):
        if r[i] < r[i-1] and r[i] < r[i+1]:
            minima.append(i)
    return minima
mins = find_local_minima(r)
# find local maxima (apoapsis)
maxs = []
for i in range(1,len(r)-1):
    if r[i] > r[i-1] and r[i] > r[i+1]:
        maxs.append(i)
# sample lengths
print("# periapsis counts:", len(mins))
print("# apoapsis counts:", len(maxs))
# sample values
sample_peri = r[mins][:10]
sample_apo = r[maxs][:10]
print("sample periapsis radii:", sample_peri)
print("sample apoapsis radii:", sample_apo)
print("sample peri times diffs (first 9):", np.diff(times[mins])[:9])
# compute mean T
diffs = np.diff(times[mins])
T_avg = np.mean(diffs)
# compute mean periapsis and apoapsis radii
r_peri_mean = np.mean(r[mins])
r_apo_mean = np.mean(r[maxs])
# a_rel
a_rel = (r_peri_mean + r_apo_mean) / 2
# compute M using Kepler's third law
G = 6.67430e-11
M = 4 * np.pi**2 * a_rel**3 / (G * T_avg**2)
m1 = m2 = M / 2
# total energy: E = - G m1 m2 / (2 a_rel)
E_total = - G * m1 * m2 / (2 * a_rel)
# print results
print("T_avg (s):", T_avg)
print("r_peri_mean (m):", r_peri_mean)
print("r_apo_mean (m):", r_apo_mean)
print("a_rel (m):", a_rel)
print("Total mass M (kg):", M)
print("m1 (kg):", m1)
print("E_total (J):", E_total)

{'answer': -5.623389804573717e+39}

import numpy as np
# Extract data from df
t = df['time'].values
star1 = df[['star1_x','star1_y','star1_z']].values
star2 = df[['star2_x','star2_y','star2_z']].values

# Compute radius of each star relative to origin
r1 = np.linalg.norm(star1, axis=1)
r2 = np.linalg.norm(star2, axis=1)
r1_mean = np.mean(r1)
r2_mean = np.mean(r2)

# Compute separation between the two stars
sep = np.linalg.norm(star2 - star1, axis=1)
a = np.mean(sep)

# Compute angle of relative position in xy-plane
rel = star2 - star1
angles = np.arctan2(rel[:,1], rel[:,0])
angles_unwrapped = np.unwrap(angles)
# Estimate angular velocity by linear fit
coeffs = np.polyfit(t, angles_unwrapped, 1)
omega = coeffs[0]  # rad/s

# Orbital period
period = 2 * np.pi / omega

# Gravitational constant
g = 6.67430e-11

# Total mass from Kepler's third law
M_total = 4 * np.pi**2 * a**3 / (g * period**2)

# Mass ratio and individual masses
M1 = M_total * (r2_mean / a)
M2 = M_total * (r1_mean / a)

# Total mechanical energy: E = -G*M1*M2/(2a)
E_total = - g * M1 * M2 / (2 * a)

# Print intermediate values and result
print("r1_mean =", r1_mean)
print("r2_mean =", r2_mean)
print("a (mean separation) =", a)
print("period =", period)
print("M_total =", M_total)
print("M1 =", M1)
print("M2 =", M2)
print("E_total =", E_total)

print(df.shape)
print(df.head())
print(df[['star1_z','star2_z']].describe())

import numpy as np
# Using previous df loaded

t = df['time'].values
x1 = df['star1_x'].values
y1 = df['star1_y'].values
z1 = df['star1_z'].values
x2 = df['star2_x'].values
y2 = df['star2_y'].values
z2 = df['star2_z'].values
# compute separation stats
sep = np.sqrt((x2-x1)**2 + (y2-y1)**2 + (z2-z1)**2)
sep_mean = np.mean(sep)
sep_std = np.std(sep)
print('sep_mean =', sep_mean)
print('sep_std =', sep_std)
print('sep_std/sep_mean =', sep_std/sep_mean)
# compute angular rate (fit) and M_total
import numpy as np
rel = np.vstack((x2-x1, y2-y1)).T
angles = np.arctan2(rel[:,1], rel[:,0])
angles_unwrapped = np.unwrap(angles)
coeffs = np.polyfit(t, angles_unwrapped, 1)
omega = coeffs[0]
period = 2*np.pi/abs(omega)
# gravitational constant
G = 6.67430e-11
# M_total
M_total = omega**2 * sep_mean**3 / G
print('omega =', omega)
print('period =', period)
print('M_total (from omega^2 a^3/G) =', M_total)
# compute dt array
dt = np.diff(t)
dt_mean = np.mean(dt)
print('dt_mean =', dt_mean, 'dt_std =', np.std(dt))
# compute velocities by central diff
# v1 at indices 1..N-2
# denom t[i+1] - t[i-1]
v1 = np.sqrt((x1[2:]-x1[:-2])**2 + (y1[2:]-y1[:-2])**2 + (z1[2:]-z1[:-2])**2) / (t[2:]-t[:-2])
v2 = np.sqrt((x2[2:]-x2[:-2])**2 + (y2[2:]-y2[:-2])**2 + (z2[2:]-z2[:-2])**2) / (t[2:]-t[:-2])
# median ratio v1/v2
ratio = np.median(v1/v2)
print('median v1/v2 =', ratio)
# mass ratio: M2/M1 = ratio
# M1 = M_total/(1+ratio), M2 = M_total*ratio/(1+ratio)
M1 = M_total/(1+ratio)
M2 = M_total*ratio/(1+ratio)
print('M1 =', M1)
print('M2 =', M2)
# E_total = -G*M1*M2/(2*sep_mean)
E_total = -G * M1 * M2 / (2 * sep_mean)
print('E_total =', E_total)

print(df[['time','star1_x','star1_y','star2_x','star2_y']].head())

# Compute min and max separation
sep = np.linalg.norm(star2-star1, axis=1)
sep_min = np.min(sep)
sep_max = np.max(sep)
print('sep_min =', sep_min)
print('sep_max =', sep_max)
print('(sep_min + sep_max)/2 =', (sep_min + sep_max)/2)
print('sep_mean =', np.mean(sep))
# Possibly average of min and max = a

import numpy as np
# positions
r1 = df[['star1_x','star1_y','star1_z']].values
r2 = df[['star2_x','star2_y','star2_z']].values
t = df['time'].values
# relative separation vector
rel = r2 - r1
sep = np.linalg.norm(rel, axis=1)
r_min = sep.min()
r_max = sep.max()
a = (r_min + r_max)/2
# angle vs time
angles = np.arctan2(rel[:,1], rel[:,0])
angles_unwrapped = np.unwrap(angles)
coeffs = np.polyfit(t, angles_unwrapped, 1)
n = abs(coeffs[0])  # mean motion
P = 2*np.pi/n
# total mass
G = 6.67430e-11
M_total = n**2 * a**3 / G
# compute v1 and v2 velocity vectors by central diff
dt = t[2:] - t[:-2]
v1_vec = (r1[2:] - r1[:-2]) / dt[:,None]
v2_vec = (r2[2:] - r2[:-2]) / dt[:,None]
# compute COM drift velocity estimate
d_v1 = v1_vec.mean(axis=0)
d_v2 = v2_vec.mean(axis=0)
print('v_cm from v1 mean =', d_v1)
print('v_cm from v2 mean =', d_v2)
v_cm = (d_v1 + d_v2)/2
print('v_cm =', v_cm)
# compute velocities relative to COM
v1_rel = v1_vec - v_cm
v2_rel = v2_vec - v_cm
# speeds
v1_rel_speed = np.linalg.norm(v1_rel, axis=1)
v2_rel_speed = np.linalg.norm(v2_rel, axis=1)
# median ratio q = m2/m1 = v1_rel_speed / v2_rel_speed
q_median = np.median(v1_rel_speed / v2_rel_speed)
# alternative average ratio
q_mean = np.mean(v1_rel_speed / v2_rel_speed)
print('q_median =', q_median)
print('q_mean =', q_mean)
# compute masses
m2_m1 = q_median  # by median
m1 = M_total/(1+m2_m1)
m2 = M_total*m2_m1/(1+m2_m1)
E_total = -G * m1 * m2 / (2 * a)
print('a =', a)
print('P =', P)
print('n =', n)
print('M_total =', M_total)
print('m1 =', m1)
print('m2 =', m2)
print('E_total =', E_total)

print('mean of star1_x:', np.mean(x1))
print('mean of star1_y:', np.mean(y1))
print('mean of star2_x:', np.mean(x2))
print('mean of star2_y:', np.mean(y2))

import numpy as np
r1 = df[['star1_x','star1_y','star1_z']].values
r2 = df[['star2_x','star2_y','star2_z']].values
t = df['time'].values
rel = r2 - r1
sep = np.linalg.norm(rel, axis=1)
r_min = sep.min()
r_max = sep.max()
a = (r_min + r_max) / 2
# compute mean motion
ingles = np.unwrap(np.arctan2(rel[:,1], rel[:,0]))
coeffs = np.polyfit(t, angles_unwrapped, 1)
n = abs(coeffs[0])
# total mass
G = 6.67430e-11
M_total = n**2 * a**3 / G
# velocity central diff
dt = t[2:] - t[:-2]
v1 = (r1[2:] - r1[:-2]) / dt[:,None]
v2 = (r2[2:] - r2[:-2]) / dt[:,None]
v1_mean = v1.mean(axis=0)
v2_mean = v2.mean(axis=0)
v_cm = 0.5*(v1_mean + v2_mean)
# relative velocity
tmp1 = v1 - v_cm
tmp2 = v2 - v_cm
speed1 = np.linalg.norm(tmp1, axis=1)
speed2 = np.linalg.norm(tmp2, axis=1)
q = np.median(speed1/speed2)
M1 = M_total / (1+q)
M2 = M_total*q/(1+q)
E_total = - G * M1 * M2 / (2 * a)
print(f"r_min={r_min}, r_max={r_max}, a={a}")
print(f"n={n}, M_total={M_total}")
print(f"q={q}, M1={M1}, M2={M2}")
print(f"E_total={E_total}")

import numpy as np
# extract arrays
t = df['time'].values
r1 = df[['star1_x','star1_y','star1_z']].values
r2 = df[['star2_x','star2_y','star2_z']].values
# compute semimajor axis a = (r_min + r_max)/2
sep = np.linalg.norm(r2 - r1, axis=1)
r_min = sep.min()
r_max = sep.max()
a = 0.5 * (r_min + r_max)
# compute mean motion n from angle unwrapped fit or period measurement
rel = r2 - r1
angles = np.arctan2(rel[:,1], rel[:,0])
angles_unwrapped = np.unwrap(angles)
# use polyfit
coeffs = np.polyfit(t, angles_unwrapped, 1)
n_fit = coeffs[0]
# but periods: absolute value
n = abs(n_fit)
P = 2 * np.pi / n
# compute q via acceleration ratio from central difference non-uniform dt
dt1 = t[1:-1] - t[:-2]
dt2 = t[2:] - t[1:-1]
dt12 = t[2:] - t[:-2]
# velocities
v1_prev = (r1[1:-1] - r1[:-2]) / dt1[:, None]
v1_next = (r1[2:] - r1[1:-1]) / dt2[:, None]
v2_prev = (r2[1:-1] - r2[:-2]) / dt1[:, None]
v2_next = (r2[2:] - r2[1:-1]) / dt2[:, None]
# accelerations
acc1 = 2 * (v1_next - v1_prev) / dt12[:, None]
acc2 = 2 * (v2_next - v2_prev) / dt12[:, None]
a1_mag = np.linalg.norm(acc1, axis=1)
a2_mag = np.linalg.norm(acc2, axis=1)
# ratio and median
ratio = a1_mag / a2_mag
q = np.median(ratio)
# gravitational constant
G = 6.67430e-11
# total mass
M_total = n**2 * a**3 / G
# individual masses
M1 = M_total / (1 + q)
M2 = M_total * q / (1 + q)
# total energy
total_energy = - G * M1 * M2 / (2 * a)
# print results
print(f"r_min = {r_min:.6e} m")
print(f"r_max = {r_max:.6e} m")
print(f"a (semi-major axis) = {a:.6e} m")
print(f"Mean motion n = {n:.6e} rad/s")
print(f"Period P = {P:.6e} s")
print(f"Mass ratio q = M2/M1 = {q:.6e}")
print(f"M_total = {M_total:.6e} kg")
print(f"M1 = {M1:.6e} kg")
print(f"M2 = {M2:.6e} kg")
print(f"Total energy = {total_energy:.6e} J")

print('t_min =', t.min(), 't_max =', t.max(), 'time span =', t.max()-t.min(), 'P=', period)

Here is the outline of how we extract the total mechanical energy of the binary system from the positional data:

1. Data and coordinate setup  
   • We have time-series positions of star 1 and star 2 in the (x,y,z) plane; z = 0 for all points so the motion is coplanar.  
   • The observations cover roughly 1.55×10¹⁰ s, which, as we will see, spans about ten orbital periods.  

2. Compute the instantaneous separation and its extrema  
   • At each time step tᵢ, form the relative vector r₂–r₁ and its magnitude  
        sᵢ = |r₂(tᵢ) – r₁(tᵢ)|.  
   • Find s_min = min sᵢ and s_max = max sᵢ.  
   • For a bound Keplerian orbit, the semimajor axis a = (s_min + s_max)/2.  
     Numerically:  
     s_min ≃ 1.875 ×10¹² m  
     s_max ≃ 7.5006 ×10¹² m  
     ⇒ a ≃ 4.6878 ×10¹² m  

3. Determine the orbital period (and mean motion)  
   • Compute the polar angle θᵢ = atan2(y₂–y₁, x₂–x₁) and unwrap it to remove 2π jumps.  
   • Fit a straight line θ_unwrapped ≃ n·t + constant by least squares.  
   • The slope gives the mean motion n ≃ 4.0312 ×10⁻⁹ rad/s, so  
        P = 2π/n ≃ 1.5586 ×10⁹ s.  

4. Total mass via Kepler’s third law  
   • Kepler’s law for two bodies:  
        a³/P² = G (M₁ + M₂)/(4π²)  
     ⇒ M_total = M₁+M₂ = 4π² a³ / (G P²)  
   • Plugging in G = 6.67430×10⁻¹¹ SI gives  
        M_total ≃ 2.5083 ×10³¹ kg.  

5. Mass ratio from relative velocities  
   • Approximate each star’s instantaneous velocity by central differences of its position.  
   • Remove the small center‐of‐mass drift to get speeds v₁, v₂ about the COM.  
   • For a two‐body orbit, v₁/v₂ = M₂/M₁ ≡ q. We find  
        q ≃ 0.32292 ⇒ M₂ ≃ 0.32292 M₁.  
   • Hence  
        M₁ = M_total/(1+q) ≃ 1.8960 ×10³¹ kg  
        M₂ = q M₁       ≃ 6.1226 ×10³⁰ kg.  

6. Total mechanical energy  
   • For a bound two‐body Kepler orbit, the sum of kinetic plus potential energy is constant and equals  
        E = − G M₁ M₂ / (2 a).  
   • Inserting our values:  
        E ≃ −(6.67430×10⁻¹¹)(1.8960×10³¹)(6.1226×10³⁰)/(2×4.6878×10¹²)  
          ≃ −8.264 ×10³⁸ J.  

Conclusion:  
The total mechanical energy of the binary system is approximately  
     E_total ≃ −8.26 ×10³⁸ joules.

Using the positional data and Kepler’s laws, we find:

– s_min ≃ 1.875 ×10¹² m  
– s_max ≃ 7.5006 ×10¹² m  
→ a = (s_min+s_max)/2 ≃ 4.6878 ×10¹² m  
– Mean motion n ≃ 4.0312 ×10⁻⁹ rad/s → P ≃ 1.5586 ×10⁹ s  
– M_total = 4π²a³/(G P²) ≃ 2.5083 ×10³¹ kg  
– Mass ratio q = M₂/M₁ ≃ 0.32292 → M₁ ≃ 1.8960 ×10³¹ kg, M₂ ≃ 6.1226 ×10³⁰ kg  
– Total energy E = –G M₁ M₂/(2a) ≃ −8.264 ×10³⁸ J

Answer: −8.264e38 J

{'answer': -8.263943373761043e+38}

import numpy as np
t = df['time'].values
r1 = df[['star1_x','star1_y','star1_z']].values
r2 = df[['star2_x','star2_y','star2_z']].values
r_rel = np.linalg.norm(r2 - r1, axis=1)
print("t range:", t[0], t[-1])
print("r_rel min, max:", r_rel.min(), r_rel.max())
print("r_rel first 10:", r_rel[:10])

df.columns

print(df.head())
print(df.info())

print(globals().keys())

r1_norm = np.linalg.norm(r1, axis=1)
 r2_norm = np.linalg.norm(r2, axis=1)
 ratio = r1_norm / r2_norm
 print(np.mean(ratio), np.std(ratio), np.min(ratio), np.max(ratio))

r1_norm = np.linalg.norm(r1, axis=1)
 r2_norm = np.linalg.norm(r2, axis=1)
 ratio = r1_norm / r2_norm
 print(np.mean(ratio), np.std(ratio), np.min(ratio), np.max(ratio))

print(r_rel[-10:])

r1_norm = np.linalg.norm(r1, axis=1)
r2_norm = np.linalg.norm(r2, axis=1)
ratio = r1_norm / r2_norm
print("ratio mean, std, min, max:", np.mean(ratio), np.std(ratio), np.min(ratio), np.max(ratio))

# Let's downsample data to speed up computations
idxs = np.linspace(0, len(t)-1, 2000).astype(int)  # 2000 samples
t_samp = t[idxs]
r1_samp = r1[idxs]
r2_samp = r2[idxs]
print(len(t_samp))

import numpy as np
import pandas as pd
from scipy.optimize import minimize_scalar

t_global = df['time'].values
r1_global = df[['star1_x', 'star1_y', 'star1_z']].values
r2_global = df[['star2_x', 'star2_y', 'star2_z']].values
# Compute velocities using central differences
# v1[i] at t[i+1]-t[i-1]
v1 = (r1_global[2:] - r1_global[:-2]) / (t_global[2:] - t_global[:-2])[:, None]
v2 = (r2_global[2:] - r2_global[:-2]) / (t_global[2:] - t_global[:-2])[:, None]
t_v = t_global[1:-1]

# Downsample
step = 20
idx_v = np.arange(0, len(t_v), step)
v1_ds = v1[idx_v]
v2_ds = v2[idx_v]
t_ds = t_v[idx_v]
print("Total velocity samples:", len(v1), "Downsampled:", len(v1_ds))

dt = np.diff(t)
print("dt stats (s): mean, std, min, max:", np.mean(dt), np.std(dt), np.min(dt), np.max(dt))

import numpy as np
t = df['time'].values
r1 = df[['star1_x','star1_y','star1_z']].values
r2 = df[['star2_x','star2_y','star2_z']].values
# relative positions
r_rel = r2 - r1
r_rel_norm = np.linalg.norm(r_rel, axis=1)
# detect periapsis and apastron
peri_idx = np.where((r_rel_norm[1:-1] < r_rel_norm[:-2]) & (r_rel_norm[1:-1] < r_rel_norm[2:]))[0] + 1
apo_idx = np.where((r_rel_norm[1:-1] > r_rel_norm[:-2]) & (r_rel_norm[1:-1] > r_rel_norm[2:]))[0] + 1
# times and radii for peri and apo
t_peri = t[peri_idx]
r_peri = r_rel_norm[peri_idx]
t_apo = t[apo_idx]
r_apo = r_rel_norm[apo_idx]
# how many peri and apo
print(len(peri_idx), len(apo_idx))
# list first few
print('first peri times, r:', list(zip(t_peri[:5], r_peri[:5])))
print('first apo times, r:', list(zip(t_apo[:5], r_apo[:5])))

import numpy as np
from scipy.optimize import minimize_scalar
# Data
t = df['time'].values
r1 = df[['star1_x','star1_y','star1_z']].values
r2 = df[['star2_x','star2_y','star2_z']].values
x1 = df['star1_x'].values
y1 = df['star1_y'].values
x2 = df['star2_x'].values
y2 = df['star2_y'].values
# compute r_rel norm
r_rel = np.linalg.norm(r2 - r1, axis=1)
# find r_min, r_max
a1 = r_rel.min()
a2 = r_rel.max()
a = 0.5*(a1 + a2)
# find peri and apo times
t_peri = t[np.where((r_rel[1:-1] < r_rel[:-2]) & (r_rel[1:-1] < r_rel[2:]))[0] + 1]
t_apo = t[np.where((r_rel[1:-1] > r_rel[:-2]) & (r_rel[1:-1] > r_rel[2:]))[0] + 1]
# compute P from peri
dts = np.diff(t_peri)
P = np.mean(dts)
# constants
G = 6.67430e-11
# compute M
a = float(a)
P = float(P)
GM = 4 * np.pi**2 * a**3 / P**2
M = GM / G
print("r_min, r_max (m):", a1, a2)
print("a (m):", a)
print("Period P (s):", P)
print("GM (m3/s2):", GM)
print("M_total (kg):", M)
# Compute q by minimizing variance of COM
# Define f(q)
# Use pixels: x1, y1, x2, y2

def f(q):
    # COM positions
    comx = (x1 + q * x2) / (1 + q)
    comy = (y1 + q * y2) / (1 + q)
    return np.var(comx) + np.var(comy)
# minimize f(q)
res = minimize_scalar(f, bounds=(0.001, 1000), method='bounded', options={'xatol':1e-12})
q = res.x
print("q (m2/m1):", q)
# compute masses
m1 = M / (1 + q)
m2 = q * m1
print("m1, m2 (kg):", m1, m2)
# reduced mass
mu = m1 * m2 / M
print("mu (kg):", mu)
# compute total energy
E_tot = - G * m1 * m2 / (2 * a)
print("Total energy E (J):", E_tot)
# Also compute E from reduced mass formula: E_tot2 = mu * (-G*M/(2a))
E_tot2 = mu * (-G * M / (2 * a))
print("Total energy check E2 (J):", E_tot2)
# evaluate & print

{'answer': -8.380600358510325e+38}

import numpy as np
# Compute sums of positions
df_sum = df[['star1_x','star2_x','star1_y','star2_y','star1_z','star2_z']].copy()
df_sum['sum_x'] = df['star1_x'] + df['star2_x']
df_sum['sum_y'] = df['star1_y'] + df['star2_y']
df_sum['sum_z'] = df['star1_z'] + df['star2_z']
# Print statistics for sums
print("sum_x mean, std:", np.mean(df_sum['sum_x']), np.std(df_sum['sum_x']))
print("sum_y mean, std:", np.mean(df_sum['sum_y']), np.std(df_sum['sum_y']))
print("sum_z mean, std:", np.mean(df_sum['sum_z']), np.std(df_sum['sum_z']))

df.head()

# Inspect df columns and few rows
print(df.columns)
print(df.head())
print(df.tail())

import numpy as np
# Compute relative separation
x1 = df['star1_x'].values
y1 = df['star1_y'].values
z1 = df['star1_z'].values
x2 = df['star2_x'].values
y2 = df['star2_y'].values
z2 = df['star2_z'].values
times = df['time'].values
# Compute separation vector and magnitude
dr = np.vstack([x2-x1, y2-y1, z2-z1])
r_sep = np.linalg.norm(dr, axis=0)
r_min = np.min(r_sep)
r_max = np.max(r_sep)
# eccentricity
e = (r_max - r_min)/(r_max + r_min)
# compute local maxima indices for r_sep
# find indices where each point greater than neighbors
# avoid ends
maxi = np.where((r_sep[1:-1] > r_sep[:-2]) & (r_sep[1:-1] > r_sep[2:]))[0] + 1
# times of maxima
t_max = times[maxi]
# Compute periods between successive maxima
periods = np.diff(t_max)
# average period
t_avg = np.mean(periods)
# Print results
print("r_min =", r_min)
print("r_max =", r_max)
print("e =", e)
print("num of maxima found =", len(maxi))
print("average period =", t_avg, " s")

import numpy as np
import scipy.stats as stats
# given arrays from previous code
# define arrays (we'll re-get them):
import pandas as pd
df_local = df.copy()
# compute a and T again:
x1 = df_local['star1_x'].values
y1 = df_local['star1_y'].values
z1 = df_local['star1_z'].values
x2 = df_local['star2_x'].values
y2 = df_local['star2_y'].values
z2 = df_local['star2_z'].values
times = df_local['time'].values
# relative separation
dr = np.vstack([x2-x1, y2-y1, z2-z1])
r_sep = np.linalg.norm(dr, axis=0)
r_min = np.min(r_sep)
r_max = np.max(r_sep)
a = 0.5 * (r_max + r_min)
# maxima of r_sep for period T
times = df_local['time'].values
# find local maxima
maxi = np.where((r_sep[1:-1] > r_sep[:-2]) & (r_sep[1:-1] > r_sep[2:]))[0] + 1
if len(maxi) < 2:
    raise ValueError("Not enough maxima to compute period")
# times of max
t_max = times[maxi]
periods = np.diff(t_max)
T = np.mean(periods)
# masses
G = 6.67430e-11
M_total = 4 * np.pi**2 * a**3 / (G * T**2)
# detrend x and y for both stars
# Fit linear regression for x1 vs time
# using numpy polyfit
coef1 = np.polyfit(times, x1, 1)
coef2 = np.polyfit(times, x2, 1)
# detrended positions
x1_d = x1 - (coef1[0]*times + coef1[1])
x2_d = x2 - (coef2[0]*times + coef2[1])
# Similarly for y
coef1_y = np.polyfit(times, y1, 1)
coef2_y = np.polyfit(times, y2, 1)
y1_d = y1 - (coef1_y[0]*times + coef1_y[1])
y2_d = y2 - (coef2_y[0]*times + coef2_y[1])
# compute stds
std_x1 = np.std(x1_d)
std_x2 = np.std(x2_d)
std_y1 = np.std(y1_d)
std_y2 = np.std(y2_d)
# compute k from x and y
k_x = std_x1 / std_x2
k_y = std_y1 / std_y2
# average k
k = (k_x + k_y)/2
# Also slope ratio of x1_d vs x2_d? But out-of-phase
# Print results
print("r_min =", r_min)
print("r_max =", r_max)
print("a =", a)
print("T =", T)
print("M_total =", M_total, "kg")
print("std_x1 =", std_x1, "std_x2 =", std_x2, "k_x =", k_x)
print("std_y1 =", std_y1, "std_y2 =", std_y2, "k_y =", k_y)
print("k =", k)

import numpy as np
G = 6.67430e-11
M_total = 2.525938722327381e31
k = 0.3229166666666641
a = 4687828369742.379
m1 = M_total/(1+k)
m2 = M_total - m1
E = - G*m1*m2/(2*a)
print(E)

{'answer': -8.38060036057417e+38}

We need to extract from the positional‐time data everything needed to get the total (K+U) of the two‐body system.  Since no masses are given explicitly, we will infer them from the orbit via Kepler’s law, then compute the binding energy.  Here is the step-by-step plan:

1. Load the DataFrame `df` (already in memory) and pull out the time array `t`, and the two 3-vectors r₁(t) and r₂(t).

2. Form the relative separation  
   r_rel(t) = r₂(t) − r₁(t)  
   and its magnitude ‖r_rel‖ as a function of time.

3. Find the orbital period P by locating successive maxima in ‖r_rel‖ (apoapsis) and taking the average Δt between peaks.

4. Estimate the semi-major axis of the relative orbit  
   a_rel ≃ (max(‖r_rel‖) + min(‖r_rel‖)) / 2.

5. Use Kepler’s third law for two bodies (μ = G(M₁+M₂)) in the reduced form  
   P² = 4π² a_rel³ / (G M_tot)  
   ⇒ M_tot = M₁ + M₂ = 4π² a_rel³ / (G P²).

6. Recover the mass ratio k = M₁/M₂ from the ratio of the individual orbital radii about the COM.  If the data are given in the COM frame, then at any instant  
   k = M₁/M₂ = ‖r₂‖ / ‖r₁‖,  
   which is constant over the whole orbit.  We can form the median of ‖r₂‖/‖r₁‖ to get a robust estimate of k.

7. From  
   M₁ = (k/(1+k)) M_tot,    M₂ = (1/(1+k)) M_tot.

8. Finally compute the total mechanical energy of the bound orbit, which for a Keplerian two‐body system is  
   E = K + U = − G M₁ M₂ / (2 a_rel).

Since in a closed Keplerian orbit the energy is constant, we only need to compute it once, at any convenient point in the data.  

Next step: implement these steps in the Python REPL, printing out   
   – the estimated P and a_rel  
   – the total mass M_tot  
   – the mass ratio k, and hence M₁, M₂  
   – the final E in joules  

Let me now drive the REPL to carry out those calculations.

import numpy as np
import pandas as pd

df.head()

print(df.columns)
print(df.head())

print(df.iloc[0].to_dict())

# Preliminary processing: compute r1, r2, r_rel, magnitudes
time = df['time'].values
r1 = df[['star1_x', 'star1_y', 'star1_z']].values
r2 = df[['star2_x', 'star2_y', 'star2_z']].values
# relative separation
r_rel = r2 - r1
# magnitudes
d_rel = np.linalg.norm(r_rel, axis=1)
r1_mag = np.linalg.norm(r1, axis=1)
r2_mag = np.linalg.norm(r2, axis=1)
# print time span and number of points and some magnitudes
print("Number of points, time span (s):", len(time), time.min(), time.max())
print("Relative distances (first few):", d_rel[:5])
print("r1 magnitudes (first few):", r1_mag[:5])
print("r2 magnitudes (first few):", r2_mag[:5])

import numpy as np
import pandas as pd
# gravitational constant
g = 6.67430e-11
# extract data
time = df['time'].values
r1 = df[['star1_x', 'star1_y', 'star1_z']].values
r2 = df[['star2_x', 'star2_y', 'star2_z']].values
# relative vector and distance
d_rel = np.linalg.norm(r2 - r1, axis=1)
# compute dt (assumes uniform sampling)
dt = time[1] - time[0]
# compute accelerations via central differences (exclude first and last)
a1 = (r1[2:] - 2*r1[1:-1] + r1[:-2]) / dt**2
a2 = (r2[2:] - 2*r2[1:-1] + r2[:-2]) / dt**2
# mid-step relative vectors and distance
r_rel_mid = r2[1:-1] - r1[1:-1]
d_rel_mid = np.linalg.norm(r_rel_mid, axis=1)
r_rel_hat_mid = r_rel_mid / d_rel_mid[:, None]
# parallel components of acceleration
a1_par = np.sum(a1 * r_rel_hat_mid, axis=1)  # equals G M2 / d^2
# for a1, acceleration in direction of r_rel gives positive when r_rel points from star1 to star2
# compute M2 estimates
M2_est = a1_par * d_rel_mid**2 / g
# a2 parallel: acceleration on star2 towards star1 i.e. along -r_rel_hat, so project a2 onto -r_rel_hat
a2_par = np.sum(a2 * (-r_rel_hat_mid), axis=1)  # G M1 / d^2
M1_est = a2_par * d_rel_mid**2 / g
# filter out bad values: keep values where |a1_par| big enough to reduce noise
mask = np.abs(a1_par) > np.percentile(np.abs(a1_par), 10)  # remove lowest 10% acceleration -> avoid near quadrature
# likewise filter where a2_par magnitude large
o_mask = np.abs(a2_par) > np.percentile(np.abs(a2_par), 10)
mask = mask & o_mask
# median mass estimates
M1_med = np.median(M1_est[mask])
M2_med = np.median(M2_est[mask])
M_tot_acc = M1_med + M2_med
# orbital separation extremes
r_min = d_rel.min()
r_max = d_rel.max()
a_rel = 0.5*(r_min + r_max)
# energy
i_M1 = M1_med
ni_M2 = M2_med
E_binding = - g * M1_med * M2_med / (2 * a_rel)
# print results
print(f"dt = {dt:.3e} s, number of valid points = {mask.sum()} of {len(a1_par)}")
print(f"Estimated M1 (from acc) = {M1_med:.3e} kg")
print(f"Estimated M2 (from acc) = {M2_med:.3e} kg")
print(f"Total mass (acc) = {M_tot_acc:.3e} kg")
print(f"r_min = {r_min:.3e} m, r_max = {r_max:.3e} m, a_rel = {a_rel:.3e} m")
print(f"Binding energy = {E_binding:.3e} J")

# find peaks in d_rel to get orbital period
peaks = np.where((d_rel[1:-1] > d_rel[:-2]) & (d_rel[1:-1] > d_rel[2:]))[0] + 1
time_peaks = time[peaks]
dt_peaks = np.diff(time_peaks)
# filter out spurious small dt - if dt < 0.5*(median), maybe remove? or keep all
# summary
print("Number of peaks found:", len(peaks))
# summary stats
dt_median = np.median(dt_peaks)
dt_mean = np.mean(dt_peaks)
dt_std = np.std(dt_peaks)
print(f"dt_peaks: mean={dt_mean:.3e} s, median={dt_median:.3e} s, std={dt_std:.3e} s")
# Let's look at first few dt_peaks
print("first 10 dt_peaks:", dt_peaks[:10])
# Now compute M_tot from Kepler:
P = dt_median
a_rel = 0.5*(d_rel.min() + d_rel.max())
M_tot_kepler = 4 * np.pi**2 * a_rel**3 / (6.67430e-11 * P**2)
print(f"Using P={P:.3e} s and a_rel={a_rel:.3e} m, M_tot_kepler = {M_tot_kepler:.3e} kg")

print("M1_med=", M1_med)
print("M2_med=", M2_med)
print("a_rel=", a_rel)
print("E_binding=", E_binding)

{'answer': -7.070610094676271e+38}

# Inspecting DataFrame columns and first few rows
print(df.columns)
print(df.head())

# inspect time range and number of rows
print(df['time'].min(), df['time'].max(), df.shape)
# print last few rows
print(df.tail())

import numpy as np
# Relative vectors
dx = df['star2_x'].values - df['star1_x'].values
dy = df['star2_y'].values - df['star1_y'].values
r = np.sqrt(dx**2 + dy**2)
a_mean = np.mean(r)
# Remove edges to avoid noise? But r constant presumably.
# Angle and time
theta = np.arctan2(dy, dx)
theta_unwrap = np.unwrap(theta)
t = df['time'].values
t_rel = t - t[0]
# Fit linear slope
coeff = np.polyfit(t_rel, theta_unwrap, 1)
slope, intercept = coeff[0], coeff[1]
P = 2*np.pi/abs(slope)
# G and total mass
G = 6.67430e-11
M_total = 4*np.pi**2 * a_mean**3 / (G * P**2)
# COM coords with equal masses
COM_x = (df['star1_x'] + df['star2_x'])/2
COM_y = (df['star1_y'] + df['star2_y'])/2
r1 = np.sqrt((df['star1_x']-COM_x)**2 + (df['star1_y']-COM_y)**2)
r2 = np.sqrt((df['star2_x']-COM_x)**2 + (df['star2_y']-COM_y)**2)
r1m, r2m = np.mean(r1), np.mean(r2)
q_est = r1m / r2m
m1 = M_total/2
m2 = M_total/2
E_tot = -G*m1*m2 / (2*a_mean)
print('a_mean (m):', a_mean)
print('P (s):', P)
print('M_total (kg):', M_total)
print('r1_mean (m):', r1m)
print('r2_mean (m):', r2m)
print('q_est:', q_est)
print('E_total (J):', E_tot)

import numpy as np
# compute r again for first few rows
dx = df['star2_x'].values - df['star1_x'].values
dy = df['star2_y'].values - df['star1_y'].values
r = np.sqrt(dx**2 + dy**2)
# Print first few r's
to_print = list(r[:10])
print('first ten r:', to_print)
print('mean r:', np.mean(r))
# Print 1st data row's dx, dy as well
dx0, dy0 = dx[0], dy[0]
print('dx0, dy0, r0:', dx0, dy0, r[0])

# print last ten r values
print(list(r[-10:]))
print('min r:', np.min(r), 'max r:', np.max(r))

import numpy as np
# compute r, dx, dy
dx = df['star2_x'].values - df['star1_x'].values
dy = df['star2_y'].values - df['star1_y'].values
r = np.sqrt(dx**2 + dy**2)
r_min = np.min(r)
r_max = np.max(r)
a_est = (r_min + r_max)/2
# compute theta and unwrap
theta = np.arctan2(dy, dx)
theta_unwrapped = np.unwrap(theta)
t0 = df['time'].values[0]
t1 = df['time'].values[-1]
time_span = t1 - t0
# compute # cycles
delta_theta = theta_unwrapped[-1] - theta_unwrapped[0]
n_cycles = delta_theta / (2*np.pi)
# period estimate
total_time = time_span
t_period = total_time / n_cycles
# other P estimation via polyfit as above
slope_poly, intercept_poly = np.polyfit(df['time'].values - t0, theta_unwrapped, 1)
P_poly = 2*np.pi/abs(slope_poly)
# G and M_total
g = 6.67430e-11
M_total_a = 4*np.pi**2 * a_est**3 / (g * t_period**2)
M_total_poly = 4*np.pi**2 * a_est**3 / (g * P_poly**2)
# masses since m1=m2 (est)
m1 = m2 = M_total_a/2
# total energy
E_total = -g * m1 * m2 /(2*a_est)
# Print results
print('r_min:', r_min)
print('r_max:', r_max)
print('a_est:', a_est)
print('n_cycles:', n_cycles)
print('time_span:', total_time)
print('t_period from cycles:', t_period)
print('P_poly:', P_poly)
print('M_total (cycles):', M_total_a)
print('M_total (polyfit):', M_total_poly)
print('E_total (J):', E_total)

# precise recomputation
import numpy as np
G = 6.67430e-11
a = (np.min(r)+np.max(r))/2
# using P_poly from above
P_poly = 2016363295.7671947
M_total_poly = 4*np.pi**2 * a**3 /(G * P_poly**2)
m1 = m2 = M_total_poly/2
E_tot = -G*m1*m2/(2*a)
print('M_total_poly:', M_total_poly)
print('E_tot:', E_tot)

# compute slope for theta vs time properly
import numpy as np
dx = df['star2_x'] - df['star1_x']
dy = df['star2_y'] - df['star1_y']
theta = np.arctan2(dy, dx)
theta_unwrapped = np.unwrap(theta.values)
time = df['time'].values - df['time'].values[0]
slope, intercept = np.polyfit(time, theta_unwrapped, 1)
print('slope:', slope)
print('P:', 2*np.pi/abs(slope))

import numpy as np
# load data from df
x1 = df['star1_x'].values
y1 = df['star1_y'].values
x2 = df['star2_x'].values
y2 = df['star2_y'].values
t = df['time'].values
# relative separation
r_vec_x = x2 - x1
r_vec_y = y2 - y1
r = np.sqrt(r_vec_x**2 + r_vec_y**2)
r_min = np.min(r)
r_max = np.max(r)
a = (r_min + r_max)/2.0
# compute theta unwrapped
theta = np.arctan2(r_vec_y, r_vec_x)
theta_unwrapped = np.unwrap(theta)
# compute time relative to first measurement
t_rel = t - t[0]
time_span = t_rel[-1]
# compute slope and poly period
slope, intercept = np.polyfit(t_rel, theta_unwrapped, 1)
P_poly = 2*np.pi/abs(slope)
# compute period from cycles
delta_theta = theta_unwrapped[-1] - theta_unwrapped[0]
n_cycles = delta_theta/(2*np.pi)
P_ncycles = abs(time_span / n_cycles)
# constants
G = 6.67430e-11
# masses and energies
def compute_M_and_E(P):
    M_tot = 4*np.pi**2 * a**3 /(G * P**2)
    m1 = m2 = M_tot/2
    E = -G * m1 * m2 / (2*a)
    return M_tot, E

M_poly, E_poly = compute_M_and_E(P_poly)
M_ncycles, E_ncycles = compute_M_and_E(P_ncycles)
print('r_min, r_max, a:', r_min, r_max, a)
sn = ' ' * 0
print(f' P_poly = {P_poly:.12e} s')
print(f' P_ncycles = {P_ncycles:.12e} s')
print('M_poly:', M_poly, 'E_poly:', E_poly)
print('M_ncycles:', M_ncycles, 'E_ncycles:', E_ncycles)

# Let's compute kinetic+potential energies per snapshot and average
import numpy as np
# Data
t = df['time'].values
x1 = df['star1_x'].values
y1 = df['star1_y'].values
dx = df['star2_x'].values - x1
dy = df['star2_y'].values - y1
r = np.sqrt(dx**2 + dy**2)
# parameters
a = (np.min(r)+np.max(r))/2.0
# compute period by cycles
theta = np.arctan2(dy, dx)
theta_unwrapped = np.unwrap(theta)
delta_theta = theta_unwrapped[-1] - theta_unwrapped[0]
n_cycles = delta_theta/(2*np.pi)
time_span = t[-1] - t[0]
P_ncycles = np.abs(time_span/n_cycles)
# masses
g = 6.67430e-11
M_total = 4*np.pi**2 * a**3 /(g * P_ncycles**2)
m1 = m2 = M_total/2.0
# compute velocities by central difference
# using vectorized array
dt = t[2:] - t[:-2]  # shape N-2
v1x = (x1[2:] - x1[:-2]) / dt
v1y = (y1[2:] - y1[:-2]) / dt
x2 = df['star2_x'].values
y2 = df['star2_y'].values
v2x = (x2[2:] - x2[:-2]) / dt
v2y = (y2[2:] - y2[:-2]) / dt
# energies
T1 = 0.5 * m1 * (v1x**2 + v1y**2)
T2 = 0.5 * m2 * (v2x**2 + v2y**2)
U = - g * m1 * m2 / r[1:-1]
E = T1 + T2 + U
# compute statistics
E_mean = np.mean(E)
E_std = np.std(E)
print('M_total:', M_total)
print('m1, m2:', m1, m2)
print('E mean:', E_mean)
print('E std dev:', E_std)
print('Percent std:', (E_std/np.abs(E_mean))*100)

import numpy as np
# Data
t = df['time'].values
x1 = df['star1_x'].values
y1 = df['star1_y'].values
x2 = df['star2_x'].values
y2 = df['star2_y'].values
# relative
dx = x2 - x1
dy = y2 - y1
r = np.sqrt(dx*dx + dy*dy)
# semi-major axis a
r_min = np.min(r)
r_max = np.max(r)
a = (r_min + r_max)/2.0
# angle vs time to get period
theta = np.arctan2(dy, dx)
theta_unwrapped = np.unwrap(theta)
t_rel = t - t[0]
slope, intercept = np.polyfit(t_rel, theta_unwrapped, 1)
P = 2*np.pi / abs(slope)
# mass via Kepler third law
g = 6.67430e-11
M_total = 4*np.pi**2 * a**3 / (g * P**2)
m1 = m2 = M_total/2.0
# total energy
e_total = - g * m1 * m2 / (2.0 * a)
# Print results
print(f"a = {a:.6e} m")
print(f"P = {P:.6e} s")
print(f"M_total = {M_total:.6e} kg")
print(f"m1 = {m1:.6e}, m2 = {m2:.6e} kg")
print(f"E_total = {e_total:.6e} J")

import numpy as np
t = df['time'].values
x1 = df['star1_x'].values;
y1 = df['star1_y'].values
x2 = df['star2_x'].values;
y2 = df['star2_y'].values
# relative
r_vec = np.sqrt((x2-x1)**2 + (y2-y1)**2)
r_min = np.min(r_vec);
r_max = np.max(r_vec)
a = (r_min + r_max)/2.0
# compute theta slope
theta = np.arctan2(y2-y1, x2-x1)
theta_unwrapped = np.unwrap(theta)
t_rel = t - t[0]
slope, _ = np.polyfit(t_rel, theta_unwrapped, 1)
P = 2*np.pi/abs(slope)
# masses and energy
G = 6.67430e-11
M_total = 4*np.pi**2 * a**3 / (G * P**2)
m1 = m2 = M_total/2.0
E_total = -G * m1 * m2/(2.0 * a)
print("r_min, r_max, a:", r_min, r_max, a)
print("P:", P)
print("M_total:", M_total)
print("E_total:", E_total)

{'answer': -9.352166084464376e+38}

import numpy as np
# Display first few rows
print(df.head())
# Compute COM assuming equal masses
com_x = (df['star1_x'] + df['star2_x']) / 2
com_y = (df['star1_y'] + df['star2_y']) / 2
com_z = (df['star1_z'] + df['star2_z']) / 2
print("COM_x sample:", com_x.head().tolist())
print("COM_y sample:", com_y.head().tolist())
print("COM_z sample:", com_z.head().tolist())
# Compute COM ranges
print("COM_x range:", com_x.max() - com_x.min())
print("COM_y range:", com_y.max() - com_y.min())
print("COM_z range:", com_z.max() - com_z.min())

# Compute time range
t0 = df['time'].min()
t1 = df['time'].max()
dt = t1 - t0
# Compute COM velocity components (approx) from computed COM_x, COM_y
com_x = (df['star1_x'] + df['star2_x']) / 2
com_y = (df['star1_y'] + df['star2_y']) / 2
vel_com_x = (com_x.max() - com_x.min()) / dt
vel_com_y = (com_y.max() - com_y.min()) / dt
# print results
print("t range (s):", dt)
print("COM v_x (m/s):", vel_com_x)
print("COM v_y (m/s):", vel_com_y)

import numpy as np
import pandas as pd
# Let's use df loaded by environment
# Compute relative separation
df_rel = df.copy()
# relative separation vector from star2 to star1
df_rel['dx'] = df_rel['star1_x'] - df_rel['star2_x']
df_rel['dy'] = df_rel['star1_y'] - df_rel['star2_y']
# z is zero
df_rel['dz'] = df_rel['star1_z'] - df_rel['star2_z']
# Compute separation
df_rel['r'] = np.sqrt(df_rel['dx']**2 + df_rel['dy']**2 + df_rel['dz']**2)
# times
t = df_rel['time'].values
r = df_rel['r'].values
# find local maxima and minima
i_max = []
i_min = []
for i in range(1, len(r)-1):
    if r[i] > r[i-1] and r[i] > r[i+1]:
        i_max.append(i)
    if r[i] < r[i-1] and r[i] < r[i+1]:
        i_min.append(i)
# times of peaks and valleys
t_peaks = t[i_max]
t_valleys = t[i_min]
r_peaks = r[i_max]
r_valleys = r[i_min]
# print counts and first few
print("Found {} peaks and {} valleys".format(len(i_max), len(i_min)))
print("First peaks times:", t_peaks[:5])
print("First valleys times:", t_valleys[:5])
# Compute periods as difference between successive peaks
T_peaks = np.diff(t_peaks)
# filter out anomalies: keep T values within some tolerance (e.g., <1.5 * median)
medT = np.median(T_peaks)
T_peaks_valid = T_peaks[np.abs(T_peaks - medT) < 0.5 * medT]
# compute average period
T = np.mean(T_peaks_valid)
# get average r_max, r_min
r_max = np.mean(r_peaks)
r_min = np.mean(r_valleys)
# Compute semi-major axis of relative orbit
a_rel = (r_max + r_min) / 2
# Report
print("Estimated period T (s):", T)
print("Estimated r_max (m):", r_max)
print("Estimated r_min (m):", r_min)
print("Estimated semi-major axis a_rel (m):", a_rel)

import numpy as np
# Extract arrays
t = df['time'].values
x1 = df['star1_x'].values;
y1 = df['star1_y'].values;
z1 = df['star1_z'].values
x2 = df['star2_x'].values;
y2 = df['star2_y'].values;
z2 = df['star2_z'].values
# Compute COM positions
com_x = (x1 + x2)/2
com_y = (y1 + y2)/2
# Fit COM drift
coeffs_x = np.polyfit(t, com_x, 1)
coeffs_y = np.polyfit(t, com_y, 1)
slope_x, intercept_x = coeffs_x;
slope_y, intercept_y = coeffs_y
# Predicted COM positions
t_com_pred = t
com_x_pred = slope_x*t + intercept_x
com_y_pred = slope_y*t + intercept_y
# Distances from COM
r1 = np.sqrt((x1-com_x_pred)**2 + (y1-com_y_pred)**2 + (z1-0)**2)
r2 = np.sqrt((x2-com_x_pred)**2 + (y2-com_y_pred)**2 + (z2-0)**2)
# relative separation
r_rel = np.sqrt((x1-x2)**2 + (y1-y2)**2 + (z1-z2)**2)
# find peaks and valleys for r_rel
i_max_rel = [i for i in range(1, len(r_rel)-1) if r_rel[i] > r_rel[i-1] and r_rel[i] > r_rel[i+1]]
i_min_rel = [i for i in range(1, len(r_rel)-1) if r_rel[i] < r_rel[i-1] and r_rel[i] < r_rel[i+1]]
# times peaks and periods
t_peaks = t[i_max_rel]
T_peaks = np.diff(t_peaks)
medT = np.median(T_peaks)
# filter valid T values
T_valid = T_peaks[np.abs(T_peaks - medT) < 0.5 * medT]
T = np.mean(T_valid)
# semimajor axis of relative orbit
a_rel = (np.mean(r_rel[i_max_rel]) + np.mean(r_rel[i_min_rel]))/2
# semimajor axes of stars around COM
a1 = (np.mean(r1[i_max_rel]) + np.mean(r1[i_min_rel]))/2
a2 = (np.mean(r2[i_max_rel]) + np.mean(r2[i_min_rel]))/2
# mass ratio q = m1/m2 = a2/a1
q = a2 / a1
# Total mass from Kepler's third law
G = 6.67430e-11
M_total = 4 * np.pi**2 * a_rel**3 / (G * T**2)
# individual masses
m1 = q / (1 + q) * M_total
m2 = 1 / (1 + q) * M_total
# total energy
E_total = -G * m1 * m2 / (2 * a_rel)
# Print results
print("T (s):", T)
print("a_rel (m):", a_rel)
print("a1 (m):", a1)
print("a2 (m):", a2)
print("a1+a2 (m):", a1+a2)
print("q (m1/m2):", q)
print("M_total (kg):", M_total)
print("m1 (kg):", m1)
print("m2 (kg):", m2)
print("E_total (J):", E_total)

{'answer': -9.19112852174083e+38}

import numpy as np

# Inspect the first few rows and shape of the data
print(df.head())
print("Data shape:", df.shape)

print(df.columns)

print([name for name in globals().keys()])

# Let's examine r1 and r2 norms for first 100 samples
import numpy as np
r1 = np.sqrt(df['star1_x']**2 + df['star1_y']**2)
r2 = np.sqrt(df['star2_x']**2 + df['star2_y']**2)
r_rel = np.sqrt((df['star1_x']-df['star2_x'])**2 + (df['star1_y']-df['star2_y'])**2)
# compute stats for first 1000
print('r1 mean/first1000:', r1.iloc[:1000].mean(), 'std:', r1.iloc[:1000].std())
print('r2 mean/first1000:', r2.iloc[:1000].mean(), 'std:', r2.iloc[:1000].std())
print('r_rel mean/first1000:', r_rel.iloc[:1000].mean(), 'std:', r_rel.iloc[:1000].std())
print('relative std (r_rel):', (r_rel.iloc[:1000].std()/r_rel.iloc[:1000].mean()))

# Compute time range and spacing
time = df['time']
print('min time:', time.min(), 'max time:', time.max())
# time span
print('time span:', time.max() - time.min())
# time difference stats
dt = time.diff().dropna()
print('dt mean, std:', dt.mean(), dt.std(), dt.min(), dt.max())

import numpy as np
from scipy.signal import argrelextrema
# constants
G = 6.67430e-11
# load time and positions
time = df['time'].values
# compute dt (mean)
dt = np.mean(np.diff(time))
# position arrays
pos1 = df[['star1_x', 'star1_y', 'star1_z']].values
pos2 = df[['star2_x', 'star2_y', 'star2_z']].values
# compute relative separation
r_rel = np.linalg.norm(pos1 - pos2, axis=1)
# find rmin, rmax, semi-major axis
rmin = r_rel.min()
rmax = r_rel.max()
a = (rmin + rmax)/2
print(f"rmin = {rmin:.6e} m, rmax = {rmax:.6e} m, a = {a:.6e} m")
# find orbital period via downsampled peaks
ds = 200
r_rel_ds = r_rel[::ds]
t = time[::ds]
# detect peaks
i_peaks = argrelextrema(r_rel_ds, np.greater, order=50)[0]
# If too few peaks, lower order
if len(i_peaks) < 3:
    i_peaks = argrelextrema(r_rel_ds, np.greater, order=25)[0]
# compute period (in seconds)
if len(i_peaks) >= 2:
    # periods in ds-sample difference units
    ds_diffs = np.diff(i_peaks)
    # convert to time: dt_ds = dt * ds
    dt_ds = dt * ds
    periods = dt_ds * ds_diffs
    T = np.mean(periods)
    print(f"Estimated orbital period T = {T:.6e} s  (~{T/ (3600*24*365.25):.3f} years), from {len(periods)} intervals")
else:
    T = None
    print(f"Could not reliably estimate orbital period, peaks found: {len(i_peaks)}")
# Estimate mass ratio m1/m2 from acceleration ratio a2/a1
# compute acceleration by central second difference
a1 = (pos1[2:] - 2*pos1[1:-1] + pos1[:-2]) / (dt**2)
a2 = (pos2[2:] - 2*pos2[1:-1] + pos2[:-2]) / (dt**2)
a1_mag = np.linalg.norm(a1, axis=1)
a2_mag = np.linalg.norm(a2, axis=1)
# ratio m1/m2 = a2/a1
ratio = a2_mag / a1_mag
# robust metrics
ratio_mean = np.mean(ratio)
ratio_median = np.median(ratio)
# filter on high a1
threshold = np.quantile(a1_mag, 0.9)
mask = a1_mag > threshold
ratio_high = ratio[mask]
ratio_mean_high = np.mean(ratio_high)
ratio_median_high = np.median(ratio_high)
print(f"m1/m2 ratio: mean = {ratio_mean:.6e}, median = {ratio_median:.6e}")
print(f"high a1 ratio: mean = {ratio_mean_high:.6e}, median = {ratio_median_high:.6e}")
# choose ratio
ratio_final = ratio_median_high
# compute total mass
M_tot = 4 * np.pi**2 * a**3 / (G * T**2)
# individual masses
m1 = (ratio_final/(1+ratio_final)) * M_tot
m2 = M_tot - m1
print(f"Total mass M = {M_tot:.6e} kg, m1 = {m1:.6e} kg, m2 = {m2:.6e} kg")
# compute total energy E = -G m1 m2/(2a)
E_tot = -G * m1 * m2 / (2 * a)
print(f"Total energy E = {E_tot:.6e} J")

# Let's compute instantaneous total mechanical energy from data for some indices
# Using m1, m2 from previous calculation
G = 6.67430e-11
m1 = m1
m2 = m2
# Prepare velocities by central difference for positions
dt = np.mean(np.diff(time))
pos1 = pos1;
pos2 = pos2
# compute velocities for indices 1 to N-2
v1 = (pos1[2:] - pos1[:-2])/(2*dt)
v2 = (pos2[2:] - pos2[:-2])/(2*dt)
# separation for indices 1 to N-2
r_rel2 = np.linalg.norm(pos1[2:] - pos2[2:], axis=1)
# compute energies for a subset of indices to reduce compute
idx = np.arange(100000, 110000) # ~100k to 110k
# but ensure idx within 0 to len(v1)-1
idx = idx[idx < len(v1)]
# compute kinetic energy and potential
KE = 0.5*m1*np.linalg.norm(v1[idx],axis=1)**2 + 0.5*m2*np.linalg.norm(v2[idx],axis=1)**2
PE = -G*m1*m2 / r_rel2[idx]
E_inst = KE + PE
print('Mean instantaneous energy (with m1,m2 from params):', np.mean(E_inst), 'std:', np.std(E_inst))
# Compare to E_tot
print('E_tot formula:', E_tot)

{'answer': -2.069779956736076e+39}

print(df.head())
print(df.info())

print('m1' in globals(), 'm2' in globals())

print(df['time'].min(), df['time'].max())
# Inspect a few sample differences
dt = df['time'].diff().dropna()
print('dt mean:', dt.mean(), 'dt std:', dt.std())

import inspect, sys
print([name for name in globals().keys() if not name.startswith('_')])

print('M1' in globals(), 'M2' in globals(), 'M_total' in globals(), 'G' in globals())

# Compute barycenter assuming m1=m2
df['com_x'] = (df['star1_x'] + df['star2_x'])/2
df['com_y'] = (df['star1_y'] + df['star2_y'])/2
df['com_z'] = (df['star1_z'] + df['star2_z'])/2
# Check first and last
print(df[['com_x','com_y','com_z']].iloc[0])
print(df[['com_x','com_y','com_z']].iloc[-1])

import numpy as np
# relative vectors
dx = df['star2_x'].values - df['star1_x'].values
dy = df['star2_y'].values - df['star1_y'].values
dz = df['star2_z'].values - df['star1_z'].values
# separation magnitude
r = np.sqrt(dx*dx + dy*dy + dz*dz)
r_min = r.min()
r_max = r.max()
a = (r_min + r_max) / 2
print('r_min =', r_min, 'r_max =', r_max, 'semi-major axis a =', a)
# period estimation via angle
t = df['time'].values
# compute angle of relative vector and unwrap
theta = np.unwrap(np.arctan2(dy, dx))
# total number of rotations
dtheta = theta[-1] - theta[0]
dt_total = t[-1] - t[0]
N_cycles = dtheta / (2 * np.pi)
P = dt_total / N_cycles
print('dtheta =', dtheta, 'dt_total =', dt_total, 'N_cycles =', N_cycles, 'P =', P)
# compute total mass
G = 6.67430e-11
M_total = 4 * np.pi**2 * a**3 / (G * P**2)
print('Total mass M_total =', M_total)
# assume equal masses M1=M2=M_total/2
# compute total energy E = -G*M1*M2/(2a) = -G*(M_total/2)^2/(2a) = -G*M_total^2/(8*a)
E = - G * M_total**2 / (8 * a)
print('Total energy E =', E)

import numpy as np, pandas as pd
# Using df loaded before
# compute dx, dy, dz
dx = df['star2_x'].values - df['star1_x'].values
dy = df['star2_y'].values - df['star1_y'].values
dz = df['star2_z'].values - df['star1_z'].values
# separation
r = np.sqrt(dx*dx + dy*dy + dz*dz)
r_min = r.min();
r_max = r.max();
a = (r_min + r_max)/2
# compute theta
theta = np.unwrap(np.arctan2(dy, dx))
dtheta = theta[-1] - theta[0]
dt_total = df['time'].values[-1] - df['time'].values[0]
N_cycles = abs(dtheta)/(2*np.pi)
P = dt_total / N_cycles
G = 6.67430e-11
M_total = 4 * np.pi**2 * a**3 / (G * P**2)
E = - G * M_total**2 / (8 * a)
print("r_min", r_min)
print("r_max", r_max)
print("a", a)
print("dtheta", dtheta)
print("N_cycles", N_cycles)
print("dt_total", dt_total)
print("P", P)
print("M_total", M_total)
print("E", E)

{'answer': -2.804183929123185e+39}

import numpy as np
from scipy.signal import find_peaks

# Extract data
t = df['time'].values
r1 = df[['star1_x','star1_y','star1_z']].values
r2 = df[['star2_x','star2_y','star2_z']].values

# relative separation
r = r2 - r1
r_mag = np.linalg.norm(r, axis=1)

# semi-major axis
a = (np.min(r_mag) + np.max(r_mag)) / 2

# find peaks of r_mag to determine period
peaks, _ = find_peaks(r_mag)
t_peaks = t[peaks]

# Compute periods between peaks
P_list = np.diff(t_peaks)
P = np.mean(P_list)

# Compute total mass via Kepler's third law
G = 6.67430e-11
M_total = 4 * np.pi**2 * a**3 / (G * P**2)

# Compute velocities via central differences
dt = np.gradient(t)
v1 = np.gradient(r1, dt, axis=0)
v2 = np.gradient(r2, dt, axis=0)
speed1 = np.linalg.norm(v1, axis=1)
speed2 = np.linalg.norm(v2, axis=1)
ratio = speed1 / speed2
ratio = ratio[np.isfinite(ratio)]
q = np.median(ratio)

# Reduced mass
mu = q/(1+q)**2 * M_total

# Total mechanical energy
E_total = - G * M_total * mu / (2 * a)

# Print results
print(a, P, M_total, q, mu, E_total)

{'answer': -2.291407544841927e+39}

import numpy as np

# Compute the separation between the two stars at each time
sep = np.sqrt((df['star1_x'] - df['star2_x'])**2 + (df['star1_y'] - df['star2_y'])**2 + (df['star1_z'] - df['star2_z'])**2)

# Find the maximum separation (apoastron)
max_sep = sep.max()

print(max_sep)

{'answer': 1274859936951.2861}

import numpy as np
# Compute separation between stars
df['sep'] = np.sqrt((df.star1_x - df.star2_x)**2 + (df.star1_y - df.star2_y)**2 + (df.star1_z - df.star2_z)**2)
# Print maximum separation (apoastron)
print(df['sep'].max())

import numpy as np
# recompute separation robustly
coords1 = df[['star1_x','star1_y','star1_z']].values
coords2 = df[['star2_x','star2_y','star2_z']].values
sep = np.linalg.norm(coords1 - coords2, axis=1)
print(sep.max())

To find the apoastron (the maximum star–star separation) I computed for each time step the Euclidean distance between the two stars:

• sep(t) = √[(x₁–x₂)² + (y₁–y₂)² + (z₁–z₂)²]  
• apoastron = maxₜ sep(t)

That yields  
apoastron ≃ 1.2748599369512861×10¹² m

I’ll now submit that as the answer.

{'answer': 1274859936951.2861}

import numpy as np
# compute distance between the two stars at each time
distances = np.sqrt((df['star1_x'] - df['star2_x'])**2 + (df['star1_y'] - df['star2_y'])**2 + (df['star1_z'] - df['star2_z'])**2)
max_dist = distances.max()
print(max_dist)

{'answer': 1274859936951.2861}

import numpy as np
# Compute the separation distance between the two stars at each time
r = np.sqrt((df['star2_x'] - df['star1_x'])**2 +
             (df['star2_y'] - df['star1_y'])**2 +
             (df['star2_z'] - df['star1_z'])**2)
# Find the maximum separation
r_max = r.max()
print(r_max)

{'answer': 1274859936951.2861}

import numpy as np
# Compute separation distances between the two stars at each recorded time
d = np.sqrt((df['star1_x'] - df['star2_x'])**2 + (df['star1_y'] - df['star2_y'])**2 + (df['star1_z'] - df['star2_z'])**2)
# Find the maximum separation (apoastron)
max_sep = d.max()
print(max_sep)

{'answer': 1274859936951.2861}

import numpy as np
# Calculate separation between the two stars at each time
r = np.sqrt((df['star1_x'] - df['star2_x'])**2 + 
             (df['star1_y'] - df['star2_y'])**2 + 
             (df['star1_z'] - df['star2_z'])**2)
# Find the maximum separation (apoastron)
max_r = r.max()
print(max_r)

{'answer': 1274859936951.2861}

import numpy as np
# Compute the separation distance between the two stars at each time step
distances = np.sqrt((df['star1_x'] - df['star2_x'])**2 + (df['star1_y'] - df['star2_y'])**2 + (df['star1_z'] - df['star2_z'])**2)
# Find the maximum distance (apoastron)
apoastron = distances.max()
print(apoastron)

{'answer': 5590630655900.134}

import numpy as np
# Compute separation distances between the two stars at each time
dists = np.sqrt((df['star1_x']-df['star2_x'])**2 + 
                (df['star1_y']-df['star2_y'])**2 + 
                (df['star1_z']-df['star2_z'])**2)
# Compute maximum distance (apoastron)
max_dist = dists.max()
print(max_dist)

{'answer': 5590630655900.134}

import numpy as np
# Compute separation at each time
distances = np.sqrt((df['star1_x'] - df['star2_x'])**2 + 
                    (df['star1_y'] - df['star2_y'])**2 + 
                    (df['star1_z'] - df['star2_z'])**2)
# Find max separation
a_max = distances.max()
print(a_max)

{'answer': 5590630655900.134}

import numpy as np
# compute separation between stars
df['sep'] = np.sqrt((df['star2_x'] - df['star1_x'])**2 + (df['star2_y'] - df['star1_y'])**2 + (df['star2_z'] - df['star1_z'])**2)
# find maximum separation (apoastron)
max_sep = df['sep'].max()
print(max_sep)

{'answer': 5590630655900.134}

import numpy as np
# Compute separation distances
d = df[['star2_x','star2_y','star2_z']].values - df[['star1_x','star1_y','star1_z']].values
dist = np.linalg.norm(d, axis=1)
max_dist = dist.max()
print(max_dist)

{'answer': 5590630655900.134}

import numpy as np
# Calculate the separation distances between the two stars
separations = np.sqrt((df['star1_x'] - df['star2_x'])**2 + (df['star1_y'] - df['star2_y'])**2 + (df['star1_z'] - df['star2_z'])**2)
# Find the maximum separation (apoastron)
apoastron = separations.max()
print(apoastron)

# Find index of maximum separation
i_max = separations.idxmax()
# Get corresponding time
t_max = df.loc[i_max, 'time']
m_min = separations.min()
i_min = separations.idxmin()
t_min = df.loc[i_min, 'time']
print("Max separation (apoastron):", separations[i_max])
print("At time:", t_max)
print("Min separation (periastron):", m_min)
print("At time:", t_min)

{'answer': 5590630655900.134}

import numpy as np
# Compute separation distances between star1 and star2
seps = np.sqrt((df['star1_x'] - df['star2_x'])**2 +
                (df['star1_y'] - df['star2_y'])**2 +
                (df['star1_z'] - df['star2_z'])**2)
# Find the maximum separation (apoastron)
max_sep = seps.max()
print(max_sep)

{'answer': 7500618110360.977}

import numpy as np
# Compute separation between the two stars at each time
df['sep'] = np.sqrt((df['star2_x'] - df['star1_x'])**2 +
                    (df['star2_y'] - df['star1_y'])**2 +
                    (df['star2_z'] - df['star1_z'])**2)
# Find maximum separation (apoastron)
max_sep = df['sep'].max()
print(max_sep)

{'answer': 7500618110360.977}

import numpy as np
# Compute distances between star1 and star2 for each time
distances = np.sqrt((df['star1_x'] - df['star2_x'])**2 + (df['star1_y'] - df['star2_y'])**2 + (df['star1_z'] - df['star2_z'])**2)
# Print the maximum distance (apoastron)
print(distances.max())

print(f"Maximum distance: {distances.max():.6e} m")

{'answer': 7500618110360.977}

import numpy as np
# Compute separation distances
r = np.sqrt((df['star2_x'] - df['star1_x'])**2 + (df['star2_y'] - df['star1_y'])**2 + (df['star2_z'] - df['star1_z'])**2)
# Find maximum separation
r_max = r.max()
print(r_max)

{'answer': 7500618110360.977}

import numpy as np
# Compute separation
 dx = df['star2_x'] - df['star1_x']
 dy = df['star2_y'] - df['star1_y']
 dz = df['star2_z'] - df['star1_z']
 sep = np.sqrt(dx**2 + dy**2 + dz**2)
 max_sep = sep.max()
 print(max_sep)

import numpy as np

dx = df['star2_x'] - df['star1_x']
dy = df['star2_y'] - df['star1_y']
dz = df['star2_z'] - df['star1_z']
separation = np.sqrt(dx**2 + dy**2 + dz**2)
max_sep = separation.max()
print(max_sep)

max_sep_value = float(separation.max())
print(max_sep_value)

{'answer': 7500618110360.977}

# Compute separation between star1 and star2 and find maximum (apoastron)
import numpy as np
# assuming df is pre-loaded
separations = np.sqrt((df['star1_x'] - df['star2_x'])**2 + (df['star1_y'] - df['star2_y'])**2 + (df['star1_z'] - df['star2_z'])**2)
max_sep = separations.max()
print(max_sep)

{'answer': 7500618110360.977}

import numpy as np
# compute separation distances between the two stars at each time
separations = np.sqrt((df['star1_x'] - df['star2_x'])**2 + (df['star1_y'] - df['star2_y'])**2 + (df['star1_z'] - df['star2_z'])**2)
# find the maximum separation (apoastron)
apo = separations.max()
print(apo)

{'answer': 10001648418449.578}

import numpy as np
# Compute separation between the two stars at each time step
df['sep'] = np.sqrt((df['star1_x'] - df['star2_x'])**2 +
                    (df['star1_y'] - df['star2_y'])**2 +
                    (df['star1_z'] - df['star2_z'])**2)
# Print maximum separation (apoastron)
print(df['sep'].max())

{'answer': 10001648418449.578}

import numpy as np
# Compute separation distances between the two stars
d = np.sqrt((df['star1_x'] - df['star2_x'])**2 + (df['star1_y'] - df['star2_y'])**2 + (df['star1_z'] - df['star2_z'])**2)
# Print the maximum distance (apoastron)
print(d.max())

{'answer': 10001648418449.578}

import numpy as np
# Compute relative separation distances between the two stars
distances = np.sqrt((df.star2_x - df.star1_x)**2 + (df.star2_y - df.star1_y)**2 + (df.star2_z - df.star1_z)**2)
# Get the maximum distance, i.e., the apoastron
max_dist = distances.max()
print(max_dist)

{'answer': 10001648418449.578}

import numpy as np
sep = np.sqrt((df['star1_x']-df['star2_x'])**2 + (df['star1_y']-df['star2_y'])**2 + (df['star1_z']-df['star2_z'])**2)
print(sep.max())

{'answer': 10001648418449.578}

import numpy as np
# Compute separation between the two stars at each time
r = np.sqrt((df['star1_x'] - df['star2_x'])**2 +
             (df['star1_y'] - df['star2_y'])**2 +
             (df['star1_z'] - df['star2_z'])**2)
# Find maximum separation (apoastron)
max_r = r.max()
print(max_r)

{'answer': 10001648418449.578}

import numpy as np
# Compute separation between the two stars at each time
sep = np.sqrt((df['star1_x'] - df['star2_x'])**2 +
              (df['star1_y'] - df['star2_y'])**2 +
              (df['star1_z'] - df['star2_z'])**2)
# Print maximum separation - apoastron
print(np.max(sep))

{'answer': 3606740160241.329}

import numpy as np
# compute relative distance between star 1 and star 2 for each time
# Assuming df is already loaded

df['distance'] = np.sqrt((df['star1_x'] - df['star2_x'])**2 + (df['star1_y'] - df['star2_y'])**2 + (df['star1_z'] - df['star2_z'])**2)

# find the maximum distance (apoastron)
max_distance = df['distance'].max()
print(max_distance)

{'answer': 3606740160241.329}

import numpy as np
# compute separation distances
distances = np.sqrt((df['star2_x'] - df['star1_x'])**2 + (df['star2_y'] - df['star1_y'])**2 + (df['star2_z'] - df['star1_z'])**2)
# find maximum
d_max = distances.max()
print(d_max)

{'answer': 3606740160241.329}

import numpy as np
# calculate separation between star1 and star2
df['sep'] = np.sqrt((df['star1_x']-df['star2_x'])**2 + (df['star1_y']-df['star2_y'])**2 + (df['star1_z']-df['star2_z'])**2)
max_sep = df['sep'].max()
print(max_sep)

{'answer': 3606740160241.329}

import numpy as np
# Compute the separation between the two stars at each time
distances = np.sqrt((df['star1_x'] - df['star2_x'])**2 +
                     (df['star1_y'] - df['star2_y'])**2 +
                     (df['star1_z'] - df['star2_z'])**2)
# Determine the apoastron (maximum separation)
apastron_distance = distances.max()
print(apastron_distance)

{'answer': 3606740160241.329}

import numpy as np
# Compute separation
diff = df[['star1_x','star1_y','star1_z']].values - df[['star2_x','star2_y','star2_z']].values
separations = np.linalg.norm(diff, axis=1)
max_sep = separations.max()
print(max_sep)

{'answer': 3606740160241.329}

import numpy as np
# Compute relative positions
r1 = df[['star1_x','star1_y','star1_z']].values
r2 = df[['star2_x','star2_y','star2_z']].values
r_rel = r2 - r1
# Compute separation
sep = np.linalg.norm(r_rel, axis=1)
# find apoastron index
i = np.argmax(sep)
# retrieving time array
t = df['time'].values
# check if we can use central difference
i, len(sep)

import numpy as np
# Extract data
t = df['time'].values
r1 = df[['star1_x','star1_y','star1_z']].values
r2 = df[['star2_x','star2_y','star2_z']].values
r = r2 - r1
# Compute separation magnitude
dist = np.linalg.norm(r, axis=1)
# Find apoastron index
i = np.argmax(dist)
# Compute relative velocity at i
N = len(t)
if 0 < i < N-1:
    dt = t[i+1] - t[i-1]
    v = (r[i+1] - r[i-1]) / dt
elif i == 0:
    dt = t[1] - t[0]
    v = (r[1] - r[0]) / dt
else:
    dt = t[-1] - t[-2]
    v = (r[-1] - r[-2]) / dt
# Compute areal velocity
areal_velocity = 0.5 * np.linalg.norm(np.cross(r[i], v))
print('Apoastron index:', i)
print('Time at apoastron:', t[i])
print('Separation at apoastron (m):', dist[i])
print('Areal velocity (m^2/s):', areal_velocity)

{'answer': 3702338275219745.0}

import numpy as np
import pandas as pd

# Ensure df is sorted by time
df_sorted = df.sort_values('time').reset_index(drop=True)

# Compute relative position vector components
dx = df_sorted['star2_x'] - df_sorted['star1_x']
dy = df_sorted['star2_y'] - df_sorted['star1_y']
dz = df_sorted['star2_z'] - df_sorted['star1_z']

# Compute squared separation and find apoastron index
df_sorted['r2'] = dx**2 + dy**2 + dz**2
idx_max = df_sorted['r2'].idxmax()

# Compute relative velocity using central differences where possible
if idx_max > 0 and idx_max < len(df_sorted) - 1:
    t_prev = df_sorted.loc[idx_max-1, 'time']
    t_next = df_sorted.loc[idx_max+1, 'time']
    r_prev = np.array([dx.iloc[idx_max-1], dy.iloc[idx_max-1], dz.iloc[idx_max-1])
    r_next = np.array([dx.iloc[idx_max+1], dy.iloc[idx_max+1], dz.iloc[idx_max+1])
    v_rel = (r_next - r_prev) / (t_next - t_prev)
else:
    if idx_max == 0:
        t_curr = df_sorted.loc[idx_max, 'time']
        r_curr = np.array([dx.iloc[idx_max], dy.iloc[idx_max], dz.iloc[idx_max]])
        t_next = df_sorted.loc[idx_max+1, 'time']
        r_next = np.array([dx.iloc[idx_max+1], dy.iloc[idx_max+1], dz.iloc[idx_max+1]])
        v_rel = (r_next - r_curr) / (t_next - t_curr)
    else:
        t_prev = df_sorted.loc[idx_max-1, 'time']
        r_prev = np.array([dx.iloc[idx_max-1], dy.iloc[idx_max-1], dz.iloc[idx_max-1]])
        t_curr = df_sorted.loc[idx_max, 'time']
        r_curr = np.array([dx.iloc[idx_max], dy.iloc[idx_max], dz.iloc[idx_max]])
        v_rel = (r_curr - r_prev) / (t_curr - t_prev)

# Compute position at apoastron
r_vec = np.array([dx.iloc[idx_max], dy.iloc[idx_max], dz.iloc[idx_max]])

# Compute rate of area swept per unit time
area_rate = 0.5 * np.linalg.norm(np.cross(r_vec, v_rel))

print(area_rate)

import numpy as np
import pandas as pd

# Ensure df is sorted by time
df_sorted = df.sort_values('time').reset_index(drop=True)

# Compute relative position vector components
dx = df_sorted['star2_x'] - df_sorted['star1_x']
dy = df_sorted['star2_y'] - df_sorted['star1_y']
dz = df_sorted['star2_z'] - df_sorted['star1_z']

# Compute squared separation and find apoastron index
r2 = dx**2 + dy**2 + dz**2
idx_max = r2.idxmax()

# Number of points
good = len(df_sorted)

# Compute relative velocity using central differences where possible
if idx_max > 0 and idx_max < good - 1:
    # central difference
    t_prev = df_sorted.loc[idx_max-1, 'time']
    t_next = df_sorted.loc[idx_max+1, 'time']
    r_prev = np.array([dx.iloc[idx_max-1], dy.iloc[idx_max-1], dz.iloc[idx_max-1]])
    r_next = np.array([dx.iloc[idx_max+1], dy.iloc[idx_max+1], dz.iloc[idx_max+1]])
    v_rel = (r_next - r_prev) / (t_next - t_prev)
elif idx_max == 0:
    # forward difference
    t_curr = df_sorted.loc[idx_max, 'time']
    t_next = df_sorted.loc[idx_max+1, 'time']
    r_curr = np.array([dx.iloc[idx_max], dy.iloc[idx_max], dz.iloc[idx_max]])
    r_next = np.array([dx.iloc[idx_max+1], dy.iloc[idx_max+1], dz.iloc[idx_max+1]])
    v_rel = (r_next - r_curr) / (t_next - t_curr)
else:
    # idx_max == last
    t_prev = df_sorted.loc[idx_max-1, 'time']
    t_curr = df_sorted.loc[idx_max, 'time']
    r_prev = np.array([dx.iloc[idx_max-1], dy.iloc[idx_max-1], dz.iloc[idx_max-1]])
    r_curr = np.array([dx.iloc[idx_max], dy.iloc[idx_max], dz.iloc[idx_max]])
    v_rel = (r_curr - r_prev) / (t_curr - t_prev)

# Position at apoastron
r_vec = np.array([dx.iloc[idx_max], dy.iloc[idx_max], dz.iloc[idx_max]])

# Compute rate of area swept per unit time
area_rate = 0.5 * np.linalg.norm(np.cross(r_vec, v_rel))

print(area_rate)

{'answer': 3702338275219745.0}

import numpy as np
# extract arrays
t = df['time'].values
r1 = df[['star1_x','star1_y','star1_z']].values
r2 = df[['star2_x','star2_y','star2_z']].values
# relative vector
r = r2 - r1
# separation magnitude for each time
sep = np.linalg.norm(r, axis=1)
# find index of maximum separation (apoastron)
i_max = np.argmax(sep)
# compute relative velocity v_rel using central difference if possible
if i_max == 0:
    v_rel = (r[1] - r[0])/(t[1] - t[0])
elif i_max == len(t)-1:
    v_rel = (r[-1] - r[-2])/(t[-1] - t[-2])
else:
    v_rel = (r[i_max+1] - r[i_max-1])/(t[i_max+1] - t[i_max-1])
# compute areal velocity at apoastron: 1/2 |r x v|
cross_vec = np.cross(r[i_max], v_rel)
dA_dt = 0.5 * np.linalg.norm(cross_vec)
# print the result
print(dA_dt)

# compute derivative using numpy.gradient
t = df['time'].values
r = df[['star1_x','star1_y','star1_z']].values
# compute v via gradient for each component
dxdt = np.gradient(r[:,0], t)
dydt = np.gradient(r[:,1], t)
dzdt = np.gradient(r[:,2], t)
v_grad = np.vstack((dxdt, dydt, dzdt)).T
# get index of max separation
i_max = np.argmax(np.linalg.norm(r, axis=1))
# compute areal velocity 1/2|r x v|
cross_vec = np.cross(r[i_max], v_grad[i_max])
dA_dt_grad = 0.5 * np.linalg.norm(cross_vec)
print(dA_dt_grad)

import numpy as np
# compute relative positions and velocities

t = df['time'].values
r1 = df[['star1_x','star1_y','star1_z']].values
r2 = df[['star2_x','star2_y','star2_z']].values
r_rel = r2 - r1
# get index of max separation
i_max = np.argmax(np.linalg.norm(r_rel, axis=1))
# compute velocities via gradient
vx = np.gradient(r_rel[:,0], t)
vy = np.gradient(r_rel[:,1], t)
vz = np.gradient(r_rel[:,2], t)
v_rel = np.vstack((vx, vy, vz)).T
# compute areal velocity
dA_dt = 0.5 * np.linalg.norm(np.cross(r_rel[i_max], v_rel[i_max]))
print(dA_dt)

{'answer': 3702338275219936.5}

import numpy as np
import pandas as pd

# Inspect the first few rows of df
def inspect_df(df):
    # Print head and tail and time differences
    print("First five rows:")
    print(df.head())
    print("\nLast five rows:")
    print(df.tail())
    print("\nTime differences:")
    dts = np.diff(df['time'].values)
    print(pd.Series(dts).describe())
    return
df = globals().get('df', None)
inspect_df(df)

import numpy as np
# get df from globals
import pandas as pd
from numpy.linalg import norm
# df is global
r = df[['star2_x','star2_y','star2_z']].values - df[['star1_x','star1_y','star1_z']].values
t = df['time'].values
# compute velocities
v = np.zeros_like(r)
# central differences and endpoints
v[1:-1] = (r[2:] - r[:-2]) / (t[2:,None] - t[:-2,None])
v[0] = (r[1] - r[0]) / (t[1] - t[0])
v[-1] = (r[-1] - r[-2]) / (t[-1] - t[-2])
# cross products
cross = np.cross(r, v)
area_rate = 0.5 * np.abs(cross[:,2])
# find apoastron index where r magnitude is max
i_max = np.argmax(np.linalg.norm(r, axis=1))
area_rate_apo = area_rate[i_max]
t_apo = t[i_max]
# print results
print("Apoastron index:", i_max)
print("Time at apoastron (s):", t_apo)
print("Max separation (m):", np.linalg.norm(r[i_max]))
print("Area rate at apoastron (m^2/s):", area_rate_apo)
print("Mean area rate (m^2/s):", area_rate.mean())
print("Std area rate (m^2/s):", area_rate.std())

{'answer': 3702338275219745.0}

import numpy as np
# Compute relative positions and times
r = df[['star2_x','star2_y','star2_z']].values - df[['star1_x','star1_y','star1_z']].values
t = df['time'].values
# Compute separation distances and find apoastron index
d = np.linalg.norm(r, axis=1)
i_apo = np.argmax(d)
# Compute velocity at apoastron using central differences if possible
if 0 < i_apo < len(t)-1:
    dt = t[i_apo+1] - t[i_apo-1]
    v = (r[i_apo+1] - r[i_apo-1]) / dt
elif i_apo == 0:
    dt = t[1] - t[0]
    v = (r[1] - r[0]) / dt
else:
    dt = t[-1] - t[-2]
    v = (r[-1] - r[-2]) / dt
# Compute areal velocity
areal_rate = 0.5 * np.linalg.norm(np.cross(r[i_apo], v))
print(areal_rate)
print("i_apo:", i_apo)
print("time_apo:", t[i_apo])

{'answer': 3702338275219745.0}

Here’s the plan to get the apoastron areal‐velocity:

1. Compute the relative position vector  
   r(t) = (x₂–x₁, y₂–y₁, z₂–z₁)  

2. Find the time index tₐ where |r(t)| is maximal (apoastron).

3. Approximate the velocity at tₐ by finite differences on r(t):  
   v(tₐ) ≃ [r(tₐ₊₁) – r(tₐ₋₁)] / [tₐ₊₁ – tₐ₋₁]  

4. Compute the instantaneous areal velocity  
   Ȧ = ½ | r(tₐ) × v(tₐ) |  

5. Report Ȧ in m²/s.

Next, I’ll load these steps into the Python REPL on our pre‐loaded DataFrame to get the numerical result.

import numpy as np
# Compute relative position vectors r(t)
r = df[['star2_x','star2_y','star2_z']].values - df[['star1_x','star1_y','star1_z']].values
# Compute their magnitudes
r_mag = np.linalg.norm(r, axis=1)
# Find apoastron index (max separation)
idx_max = int(np.argmax(r_mag))
# Extract times array
t = df['time'].values
# Compute velocity via finite differences
if idx_max == 0:
    dr = r[1] - r[0]
    dt = t[1] - t[0]
elif idx_max == len(r) - 1:
    dr = r[-1] - r[-2]
    dt = t[-1] - t[-2]
else:
    dr = r[idx_max + 1] - r[idx_max - 1]
    dt = t[idx_max + 1] - t[idx_max - 1]
v = dr / dt
# Compute areal velocity A_dot = 0.5 * |r x v|
cross = np.cross(r[idx_max], v)
A_dot = 0.5 * np.linalg.norm(cross)
print(A_dot)

{'answer': 3702338275219745.0}

import numpy as np
# Compute relative position vector r for each time step
r = df[['star2_x','star2_y','star2_z']].values - df[['star1_x','star1_y','star1_z']].values
# Compute separation distances
dist = np.linalg.norm(r, axis=1)
# Find index of maximum separation (apoastron)
i_max = np.argmax(dist)
# Extract times array
t = df['time'].values
# Compute velocity vector at apoastron using finite difference
n = len(r)
if 0 < i_max < n-1:
    v = (r[i_max+1] - r[i_max-1]) / (t[i_max+1] - t[i_max-1])
elif i_max == 0:
    v = (r[1] - r[0]) / (t[1] - t[0])
else:
    v = (r[-1] - r[-2]) / (t[-1] - t[-2])
# Position at apoastron
r_vec = r[i_max]
# Compute areal velocity: 0.5 * |r x v|
cross_vec = np.cross(r_vec, v)
A_dot = 0.5 * np.linalg.norm(cross_vec)
# Print the result
print(A_dot)

{'answer': 4.255275321042627e+16}

import numpy as np
# Load df from environment
# Compute relative position vectors and times
# We'll assume df is available
r = df[['star2_x', 'star2_y', 'star2_z']].values - df[['star1_x', 'star1_y', 'star1_z']].values
time = df['time'].values
# Compute separation magnitudes
dist = np.linalg.norm(r, axis=1)
# Find index of maximum separation (apoastron)
i_max = np.argmax(dist)
# Print i_max and time at apoastron
print("i_max:", i_max)
print("time at apoastron:", time[i_max])

import numpy as np
r = df[['star2_x','star2_y','star2_z']].values - df[['star1_x','star1_y','star1_z']].values
time = df['time'].values
i_max = np.argmax(np.linalg.norm(r, axis=1))
# Compute relative velocity via central difference
if i_max > 0 and i_max < len(time)-1:
    dt = time[i_max+1] - time[i_max-1]
    v = (r[i_max+1] - r[i_max-1]) / dt
else:
    v = None
# Compute areal velocity
cross = np.cross(r[i_max], v)
area_rate = 0.5 * np.linalg.norm(cross)
print(area_rate)

# Check time spacing
dt_prev = time[i_max] - time[i_max-1]
dt_next = time[i_max+1] - time[i_max]
print("dt_prev, dt_next:", dt_prev, dt_next)
# Check few dt uniformity
dts = np.diff(time)
print("Mean dt:", np.mean(dts), "Std dt:", np.std(dts))

{'answer': 4.255275321042627e+16}

import numpy as np
# Extract necessary arrays from df
star1 = df[['star1_x','star1_y','star1_z']].values
star2 = df[['star2_x','star2_y','star2_z']].values
times = df['time'].values
# Compute relative position vectors
r = star2 - star1# shape (n,3)
# Distances
distances = np.linalg.norm(r, axis=1)
# Index of apoastron (max distance)
idx = np.argmax(distances)
# Compute velocity at idx using central difference if possible
if 0 < idx < len(times)-1:
    dt = times[idx+1] - times[idx-1]
    v = (r[idx+1] - r[idx-1]) / dt
else:
    # use numpy gradient for edge
    v = np.gradient(r, times, axis=0)[idx]
# Compute areal sweep rate
areal_rate = 0.5 * np.linalg.norm(np.cross(r[idx], v))
print(areal_rate)

{'answer': 4.255275321042627e+16}

import numpy as np
# Compute relative position vectors
r = df[['star2_x','star2_y','star2_z']].values - df[['star1_x','star1_y','star1_z']].values
# Time array
 t = df['time'].values
# Compute separation magnitudes
sep = np.linalg.norm(r, axis=1)
# Find index of apoastron (maximum separation)
i_max = np.argmax(sep)
# Compute relative velocity at apoastron using central difference
if 0 < i_max < len(t) - 1:
    v_i = (r[i_max + 1] - r[i_max - 1]) / (t[i_max + 1] - t[i_max - 1])
elif i_max == 0:
    v_i = (r[1] - r[0]) / (t[1] - t[0])
else:
    v_i = (r[-1] - r[-2]) / (t[-1] - t[-2])
# Position vector at apoastron
r_i = r[i_max]
# Areal velocity: 0.5 * |r x v|
areal_velocity = 0.5 * np.linalg.norm(np.cross(r_i, v_i))
print(i_max, t[i_max], areal_velocity)

import numpy as np
# Compute relative position vectors
r = df[['star2_x','star2_y','star2_z']].values - df[['star1_x','star1_y','star1_z']].values
# Time array
t = df['time'].values
# Compute separation magnitudes
sep = np.linalg.norm(r, axis=1)
# Find index of apoastron (maximum separation)
i_max = np.argmax(sep)
# Compute relative velocity at apoastron using central difference
if 0 < i_max < len(t) - 1:
    dt = t[i_max + 1] - t[i_max - 1]
    v_i = (r[i_max + 1] - r[i_max - 1]) / dt
elif i_max == 0:
    v_i = (r[1] - r[0]) / (t[1] - t[0])
else:
    v_i = (r[-1] - r[-2]) / (t[-1] - t[-2])
# Position vector at apoastron
r_i = r[i_max]
# Areal velocity: 0.5 * |r x v|
areal_velocity = 0.5 * np.linalg.norm(np.cross(r_i, v_i))
print(i_max, t[i_max], areal_velocity)

{'answer': 4.255275321042627e+16}

import numpy as np
# Compute relative position vector between star2 and star1
r_vec = np.vstack((df['star2_x'] - df['star1_x'],
                   df['star2_y'] - df['star1_y'],
                   df['star2_z'] - df['star1_z'])).T  # shape (N,3)

# Compute separation magnitude at each time step
r_mag = np.linalg.norm(r_vec, axis=1)

# Find index of maximum separation (apoastron)
i_max = np.argmax(r_mag)

# Print index and corresponding time for reference
print("Index of apoastron (max separation):", i_max)
print("Time at apoastron [s]:", df['time'].values[i_max])

# Compute velocity at i_max via central difference
t = df['time'].values
N = len(t)

# Decide difference method
if 0 < i_max < N-1:
    # central difference
dt = t[i_max+1] - t[i_max-1]
v_i = (r_vec[i_max+1] - r_vec[i_max-1]) / dt
else:
    # forward/backward difference
    if i_max == 0:
        dt = t[1] - t[0]
        v_i = (r_vec[1] - r_vec[0]) / dt
    else:
        dt = t[-1] - t[-2]
        v_i = (r_vec[-1] - r_vec[-2]) / dt

# Relative position at apoastron
a_r = r_vec[i_max]

# Compute area sweep rate: (1/2) * |r x v|
cross = np.cross(a_r, v_i)
area_rate = 0.5 * np.linalg.norm(cross)

print("Area sweep rate at apoastron [m^2/s]:", area_rate)

import numpy as np
# Compute relative position vector between star2 and star1
r_vec = np.vstack((df['star2_x'] - df['star1_x'],
                   df['star2_y'] - df['star1_y'],
                   df['star2_z'] - df['star1_z'])).T  # shape (N,3)

# Compute separation magnitude at each time step
r_mag = np.linalg.norm(r_vec, axis=1)
# Find index of maximum separation (apoastron)
i_max = np.argmax(r_mag)
time_apo = df['time'].values[i_max]

# Compute velocity at i_max via finite difference
t = df['time'].values
N = len(t)
if 0 < i_max < N-1:
    # central difference
    dt = t[i_max+1] - t[i_max-1]
    v_i = (r_vec[i_max+1] - r_vec[i_max-1]) / dt
elif i_max == 0:
    dt = t[1] - t[0]
    v_i = (r_vec[1] - r_vec[0]) / dt
else:
    dt = t[-1] - t[-2]
    v_i = (r_vec[-1] - r_vec[-2]) / dt

# Relative position at apoastron
a_r = r_vec[i_max]
# Compute areal velocity: (1/2) * |r x v|
cross = np.cross(a_r, v_i)
area_rate = 0.5 * np.linalg.norm(cross)

print("Index of apoastron:", i_max)
print("Time at apoastron [s]:", time_apo)
print("Area sweep rate at apoastron [m^2/s]:", area_rate)

{'answer': 4.255275321042627e+16}

import numpy as np
# compute relative positions
r = df[['star2_x','star2_y','star2_z']].values - df[['star1_x','star1_y','star1_z']].values
# times
t = df['time'].values
# distances
dist = np.linalg.norm(r, axis=1)
# index of max r (apoastron)
i_max = np.argmax(dist)
# relative velocity at apoastron via finite difference
if i_max == 0:
    v = (r[1] - r[0]) / (t[1] - t[0])
elif i_max == len(r) - 1:
    v = (r[-1] - r[-2]) / (t[-1] - t[-2])
else:
    v = (r[i_max + 1] - r[i_max - 1]) / (t[i_max + 1] - t[i_max - 1])
# position at apoastron
r_i = r[i_max]
# cross product and area rate
dA_dt = 0.5 * np.linalg.norm(np.cross(r_i, v))
print(dA_dt)

# Recompute and print details for verification
import numpy as np
# relative positions
r = df[['star2_x','star2_y','star2_z']].values - df[['star1_x','star1_y','star1_z']].values
# times
t = df['time'].values
# distances
orm_r = np.linalg.norm(r, axis=1)
i_max = np.argmax(norm_r)
# compute velocity via central difference (or forward/back)
if i_max == 0:
    v = (r[1] - r[0]) / (t[1] - t[0])
elif i_max == len(r) - 1:
    v = (r[-1] - r[-2]) / (t[-1] - t[-2])
else:
    v = (r[i_max + 1] - r[i_max - 1]) / (t[i_max + 1] - t[i_max - 1])
r_i = r[i_max]
cross = np.cross(r_i, v)
dA_dt = 0.5 * np.linalg.norm(cross)
print("i_max:", i_max)
print("t[i_max]:", t[i_max])
print("r_i norm (distance):", norm_r[i_max])
print("v (velocity vector):", v)
print("|v|:", np.linalg.norm(v))
print("cross norm:", np.linalg.norm(cross))
print("dA/dt:", dA_dt)

import numpy as np
# Compute the area sweep rate at apoastron
r = df[['star2_x','star2_y','star2_z']].values - df[['star1_x','star1_y','star1_z']].values
# times
t = df['time'].values
# compute distances
norms = np.linalg.norm(r, axis=1)
# find apoastron index
i_max = np.argmax(norms)
# compute velocity via finite difference
if i_max == 0:
    v = (r[1] - r[0]) / (t[1] - t[0])
elif i_max == len(r) - 1:
    v = (r[-1] - r[-2]) / (t[-1] - t[-2])
else:
    v = (r[i_max+1] - r[i_max-1]) / (t[i_max+1] - t[i_max-1])
# relative position at apoastron
r_i = r[i_max]
# cross product
cross = np.cross(r_i, v)
# area rate = 1/2 |cross|
dA_dt = 0.5 * np.linalg.norm(cross)

# print results
print("apoastron index:", i_max)
print("time at apoastron (s):", t[i_max])
print("distance at apoastron (m):", norms[i_max])
print("velocity magnitude at apoastron (m/s):", np.linalg.norm(v))
print("cross product magnitude:", np.linalg.norm(cross))
print("area sweep rate (m^2/s):", dA_dt)

# Check time spacing around apoastron
ti = t[i_max]
dt1 = ti - t[i_max-1]
dt2 = t[i_max+1] - ti
# Check default dt for first few points
dts = np.diff(t[:10])
print('dt1:', dt1)
print('dt2:', dt2)
print('first few dt:', dts)
print('mean dt:', np.mean(np.diff(t)))

Here’s the step-by-step outline of what we did:

1. Extract the 3D position vectors of star 2 relative to star 1 at each time:
   r(t) = (x₂–x₁, y₂–y₁, z₂–z₁).

2. Compute the separation ∥r(t)∥ at each time and find the index i_max where this is largest (apoastron).

3. Estimate the relative velocity v at that index by a centered finite difference:
   v ≈ [r(t₊) – r(t₋)] / [t₊ – t₋].

4. Compute the instantaneous areal sweep rate via Kepler’s second‐law analogue:
   dA/dt = ½‖r × v‖

5. Plug in the values at apoastron to get  
   dA/dt ≃ 4.2553×10¹⁶ m²/s.

I will now submit this result.

{'answer': 4.255275321042627e+16}

import numpy as np
# Step 1: relative position
r = df[['star2_x','star2_y','star2_z']].values - df[['star1_x','star1_y','star1_z']].values
# time array
t = df['time'].values
# Step 2: compute magnitude & find apoastron index
mag = np.linalg.norm(r, axis=1)
i0 = np.argmax(mag)
# Step 3: approximate velocity via gradient
# dt for each point
dt = np.gradient(t)
dr = np.gradient(r, axis=0)
v = dr / dt[:, None]
# Step 4: cross product and areal rate
cross = np.cross(r[i0], v[i0])
areal_rate = 0.5 * np.linalg.norm(cross)
# Print result
print(areal_rate)

{'answer': 3.555926193843243e+16}

import numpy as np
# Compute relative position vectors
r = df[['star2_x','star2_y','star2_z']].values - df[['star1_x','star1_y','star1_z']].values
# Compute relative distances
distances = np.linalg.norm(r,axis=1)
# Identify apoastron index
idx_apo = np.argmax(distances)
# Compute time array
time = df['time'].values
# Compute relative velocities by central difference
dv = np.gradient(r, time, axis=0)
# Extract r and v at apoastron
r_apo = r[idx_apo]
v_apo = dv[idx_apo]
# Compute specific angular momentum vector
h = np.cross(r_apo, v_apo)
# Rate of area sweep
dA_dt = np.linalg.norm(h)/2
print(dA_dt)

{'answer': 3.5559261938432424e+16}

import numpy as np
# compute r vectors
t = df['time'].values
r = df[['star2_x','star2_y','star2_z']].values - df[['star1_x','star1_y','star1_z']].values
# separation magnitudes
sep = np.linalg.norm(r,axis=1)
i = np.argmax(sep)
# compute relative velocity by central difference
t_len = len(t)
if i>0 and i<t_len-1:
    v = (r[i+1] - r[i-1])/(t[i+1] - t[i-1])
elif i == 0:
    v = (r[1] - r[0])/(t[1] - t[0])
else:
    v = (r[-1] - r[-2])/(t[-1] - t[-2])
# cross product and areal velocity
cross = np.cross(r[i], v)
areal_rate = 0.5 * np.linalg.norm(cross)
print("Apoastron index:", i)
print("Apoastron time:", t[i])
print("Areal rate at apoastron (m^2/s):", areal_rate)

# Let's sample indices
import numpy as np
N = len(t)
# pick 200 random interior points
time_indices = np.random.choice(np.arange(1, N-1), size=200, replace=False)
r_all = r
rates = []
for idx in time_indices:
    dt = t[idx+1] - t[idx-1]
    v_tmp = (r_all[idx+1] - r_all[idx-1]) / dt
    cross_tmp = np.cross(r_all[idx], v_tmp)
    rates.append(0.5 * np.linalg.norm(cross_tmp))
rates = np.array(rates)
print('mean:', rates.mean())
print('median:', np.median(rates))
print('std:', rates.std())
print('min:', rates.min(), 'max:', rates.max())

# compute alternative v via v2-v1
r1 = df[['star1_x','star1_y','star1_z']].values
r2 = df[['star2_x','star2_y','star2_z']].values
# using index i above (apoastron index)
i, areal_rate_orig = int(i), areal_rate
# But we need to fetch i from previous? We can recompute i.
sep = np.linalg.norm(r,axis=1)
i = np.argmax(sep)
# central difference for r1 and r2
if i>0 and i< len(t)-1:
    v1 = (r1[i+1] - r1[i-1])/(t[i+1] - t[i-1])
    v2 = (r2[i+1] - r2[i-1])/(t[i+1] - t[i-1])
else:
    # fallback
    v1 = (r1[1] - r1[0])/(t[1] - t[0])
    v2 = (r2[1] - r2[0])/(t[1] - t[0])
v_rel = v2 - v1
cross_alt = np.cross(r[i], v_rel)
areal_rate_alt = 0.5 * np.linalg.norm(cross_alt)
print("Original:", areal_rate)
print("Alternative:", areal_rate_alt)
print("Difference:", abs(areal_rate - areal_rate_alt))

{'answer': 3.555926193843243e+16}